diff --git a/doc/elementlibmanual/elementlibmanual.tex b/doc/elementlibmanual/elementlibmanual.tex index 54233102d..696f9e144 100644 --- a/doc/elementlibmanual/elementlibmanual.tex +++ b/doc/elementlibmanual/elementlibmanual.tex @@ -437,7 +437,7 @@ \subsubsection{Lattice3dBoundary element} \elementDescription{Reference}{\cite{AthWheGra18}} \end{elementsummary} -\subsubsection{Latticelink3d} +\subsubsection{Latticelink3d element} This element represents a two-node 3d link element connecting 3d beam and 3d lattice elements. Each node has six degrees of freedom. The input parameters for this element are shown in Table~\ref{latticelink3dsummary}. @@ -451,7 +451,7 @@ \subsubsection{Latticelink3d} \elementDescription{Reference}{\cite{GraAnt19}} \end{elementsummary} -\subsubsection{Latticelink3dboundary} +\subsubsection{Latticelink3dboundary element} Represents three-node 3d boundary link element connecting 3d beam and 3d lattice elements. The first two nodes have the same meaning as for latticelink3d. The third node is used to control the loading of the periodic cell using three normal (xx, yy and zz) and three shear strain (yz, xz, xy) components. The specific input parameters for this element in addition of those for latticelink3d are shown in Table~\ref{latticelink3dboundarysummary}. @@ -459,6 +459,33 @@ \subsubsection{Latticelink3dboundary} \elementParam{\param{location}: } \elementDescription{Unknowns}{Six dofs ($u$-displacement, $v$-displacement, $w$-displacement, $u$-rotation, $v$-rotation and $w$-rotation) are required in each node.} \elementDescription{Reference}{\cite{GraAnt19}} +\end{elementsummary} + + \subsubsection{Latticeframe3d element} + Latticeframe3d represents a two-node linear 3d frame element. Each node has six degrees of freedom, namely three translations and three rotations. The element is based on the link between rigid body spring and Timoshenko frame models. The main idea is to model the elastic and inelastic response of a connection of two nodes by a set of springs located at the contact facet of two rigid bodies. Displacement jumps are computed at the contact facet, which are smeared out over the element length in the form of strains. The input parameters for this element are shown in Table ~\ref{latticeframe3dsummary}. For the computation of the nodal forces, the initial configuration is used. +%% \begin{figure}[htb] +%% \centering +%% \includegraphics[width=0.3\textwidth]{./latticeframe.png} +%% \includegraphics[width=0.3\textwidth]{./latticeframe1.png} +%% \caption{latticeframe3d element. Node numbering, DOF numbering, cross-section vertices and local coordinate system at integration point.} +%% \label{latticeframe3dfig} +%% \end{figure} + + \begin{elementsummary}{latticeframe3d}{3d lattice frame element}{{ \field{zaxis}{ra} }{rn}}{latticeframe3d element summary}{latticeframe3dsummary} +\elementParam{\param{zaxis}: z-axis of local coordinate system } +\elementDescription{CS properties}{Area, inertia moment along y and z axis (\param{iy} and \param{iz} parameters), torsion inertia moment (\param{ik} parameter) and either cross section area shear correction factor (\param{beamshearcoeff} parameter) or equivalent shear areas (\param{shearareay} and \param{shearareaz} parameters) are required. These cross section properties are assumed to be defined in local coordinate system of element.} +\elementDescription{Unknowns}{Six dofs ($u$-displacement, $v$-displacement, $w$-displacement, $u$-rotation, $v$-rotation and $w$-rotation) are required in each node.} +\elementDescription{Reference}{\cite{Toi91, Toi93}} +\end{elementsummary} + + \subsubsection{Latticeframe3dnl element} + Latticeframe3dnl represents a two-node nonlinear 3d frame element. Each node has six degrees of freedom, namely three translations and three rotations. The element is based on the link between rigid body spring and Timoshenko frame models. The main idea is to model the elastic and inelastic response of a connection of two nodes by a set of springs located at the contact facet of two rigid bodies. Displacement jumps are computed at the contact facet considering the incrementally deformed geometry, which are smeared out over the element length in the form of strains. The input parameters for this element are shown in Table ~\ref{latticeframe3dnlsummary}. For the computation of the nodal forces, the incrementally deformed configuration is used. At the start of each step the local co-ordinate system is updated with the displacements of the previous step. + + \begin{elementsummary}{latticeframe3dnl}{3d lattice frame element}{{ \field{zaxis}{ra} }{rn}}{latticeframe3d element summary}{latticeframe3dnlsummary} +\elementParam{\param{zaxis}: z-axis of local coordinate system } +\elementDescription{CS properties}{Area, inertia moment along y and z axis (\param{iy} and \param{iz} parameters), torsion inertia moment (\param{ik} parameter) and either cross section area shear correction factor (\param{beamshearcoeff} parameter) or equivalent shear areas (\param{shearareay} and \param{shearareaz} parameters) are required. These cross section properties are assumed to be defined in local coordinate system of element.} +\elementDescription{Unknowns}{Six dofs ($u$-displacement, $v$-displacement, $w$-displacement, $u$-rotation, $v$-rotation and $w$-rotation) are required in each node.} +\elementDescription{Reference}{\cite{Toi91, Toi93, AbdGra23}} \end{elementsummary} %----------------------------------------------------------------------------------------------- @@ -2240,17 +2267,15 @@ \subsubsection{Tet1\_3D\_SUPG element} %\printbibliography \begin{thebibliography}{99} - - \bibitem{RobertCook1989} R. D. Cook and D. S. Malkus and M. E. Plesha, ``Concepts and Applications of Finite Element Analysis'', Third Edition, isbn: 0-471-84788-7, 1989. \bibitem{BittnarSejnoha1996} Z. Bittnar and J. Sejnoha, ``Numerical Methods in Structural Mechanics'',Thomas Telford,isbn:978-0784401705, 1996. \bibitem{RagnarLarsson2011} R. Larsson and J. Mediavilla and M. Fagerström, ``Dynamic fracture modeling in shell structures based on XFEM'', International Journal for Numerical Methods in Engineering, vol. 86, no. 4-5, 499--527, 2011. \bibitem{GraJir10} P. Grassl and M. Jir\'{a}sek, ``Meso-scale approach to modelling the fracture process zone of concrete subjected to uniaxial tension'', International Journal of Solids and Structures, vol. 47, iss. 7-8, pp. 957-968, 2010.. \bibitem{GraBol16} P. Grassl, J. Bolander, ``Three-Dimensional Network Model for Coupling of Fracture and Mass Transport in Quasi-Brittle Geomaterials'', Materials, 9, 782, 2016 \bibitem{AthWheGra18} I. Athanasiadis, S. Wheeler and P. Grassl. ``Hydro-mechanical network modelling of particulate composites'', International Journal of Solids and Structures, vol. 130-131, pp. 49-60, 2018. -\bibitem{GraAnt19} P. Grassl and A. Antonelli. ``3D network modelling of fracture processes in fibre-reinforced geomaterials'', International Journal of Solids and Structures, vol. 156-157, Pages 234-242, 2019. - - +\bibitem{GraAnt19} P. Grassl and A. Antonelli. ``3D network modelling of fracture processes in fibre-reinforced geomaterials'', International Journal of Solids and Structures, vol. 156-157, pp. 234-242, 2019. +\bibitem{Toi91} Y. Toi. ``Shifted integration technique in one‐dimensional plastic collapse analysis using linear and cubic finite elements'', International Journal for Numerical Methods in Engineering 31, no. 8, pp. 1537-1552, 1991. +\bibitem{AbdGra23} G. Abdelrhim and P. Grassl. ``On a geometric nonlinear Timoshenko frame element with material nonlinearity'', Annual Conference of the UK Association for Computational Mechanics (UKACM), The University of Warwick, UK, 19-21 April 2023. \end{thebibliography} \end{document} diff --git a/doc/matlibmanual/matlibmanual.tex b/doc/matlibmanual/matlibmanual.tex index bafc5e089..14aeb36d3 100644 --- a/doc/matlibmanual/matlibmanual.tex +++ b/doc/matlibmanual/matlibmanual.tex @@ -6079,6 +6079,74 @@ \subsubsection{Viscoelastic plasticity damage lattice model} \label{latticeplasticdamageviscoelastic_table} \end{table} +\subsubsection{Plasticity lattice model for steel} + +This is a Plasticity material used together with lattice elements.It has bean developed to describe the failure process failure of steel frame member. +The stress-strain law has the form +\begin{equation} +\boldsymbol{\sigma} = \mathbf{D}_{\rm e} \left(\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}_{\rm p} \right) +\end{equation} +where $\boldsymbol{\sigma}$ is a vector of tractions and rotational components, $\mathbf{D}_{\rm e}$ is the elastic stiffness matrix, $\boldsymbol{\varepsilon}$ is a vector of strains obtained from displacement jumps smeared over the element length and rotational components and $\boldsymbol{\varepsilon}_{\rm p}$ are the plastic strains. The model parameters are summarised in Table \ref{latticeframesteelplastic_table}. +\begin{table}[!htb] + \centering +\begin{tabular}{|l|l|} +\hline +Description & Plasticity lattice model for steel \\ +\hline +Record Format & \descitem{latticeframesteelplastic} $\elemparam{}_{(in)}$ +$\elemparam{d}_{(rn)}$ $\optelemparam{talpha}_{(rn)}$ $\elemparam{e}_{(rn)}$ $\optelemparam{n}_{(rn)}$ $\optelemparam{nx0}_{(rn)}$ $\optelemparam{mx0}_{(rn)}$\\ &$\optelemparam{my0}_{(rn)}$ $\optelemparam{mz0}_{(rn)}$ $\elemparam{sub}_{(in)}$ $\elemparam{iter}_{(in)}$ $\elemparam{tol}_{(rn)}$\\ +Parameters &- \param{} material number\\ +&- $\param{d}$ material density\\ +&- $\optparam{talpha}$ Thermal expansion coefficient. Default is 0.\\ +&- $\optparam{e}$ Young's modulus of the lattice material.\\ +&- $\param{n}$ Poisson's ratio of the material that the beam element is made of.\\ +&- $\optparam{nx0}$ ultimate capacity under pure axialloads.\\ +&- $\optparam{mx0}$ ultimate capacity under pure moment about $x$.\\ +&- $\optparam{my0}$ ultimate capacity under pure moment about $y$.\\ +&- $\optparam{mz0}$ ultimate capacity under pure moment about $z$.\\ +&- $\optparam{sub}$ maximum number of subincrementations.\\ +&- $\optparam{iter}$ maximum number of newton iterations.\\ +&- $\optparam{tol}$ tolerance for newton method.\\ +Supported modes & 3dlattice\\ + +\hline +\end{tabular} + +\caption{Plasticity steel model for lattice elements -- summary.} +\label{latticeframesteelplastic_table} +%\end{center} + \end{table} +The nonlinear response of the beam element is given directly as an relationship between internal force and strain. The basic equations include an additive decomposition of total strain into elastic part and plastic part + +\begin{equation} +\boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}_{\rm e}-\boldsymbol{\varepsilon}_{\rm p} +\end{equation} + +Trail stress +\begin{equation} + \sigma_{n+1}^{tr} =D^e (\varepsilon_{n+1} - \varepsilon_n^p) +\end{equation} + +The yield function +\begin{equation} +(\frac{N_x}{N_0})^2+ (\frac{M_x}{M_{x0}} )^2+ (\frac{M_y}{M_{y0}} )^2+ (\frac{M_z}{M_{z0}} )^2- 1 = 0 +\end{equation} + +where $M$, $M y$ , $N$ and $M$ , are the two components of bending moment, an axial force and +a torsional moment, respectively. The subscript 0 indicates a fully plastic value under the +condition that each component of resultant forces acts independently on a cross-section of +a member. +\\ +loading-unloading conditions + +loading-unloading conditions +\begin{equation} +f(\bar{\vsig},\kappa)\le 0 \qquad \dot{\lambda}\geq 0 \qquad \dot{\lambda}f(\bar{\vsig},\kappa)=0, +\end{equation} +evolution law for plastic strain +\begin{equation} +\dot{\veps}_{\rm{p}} = \dot{\lambda} \frac{\partial f}{\partial \bar{\vsig}}, +\end{equation} \subsection{Material models for steel relaxation} @@ -6191,7 +6259,6 @@ \subsubsection{Model for relaxation of prestressing steel - SteelRelaxMat} \end{table} - \subsection{User-defined material models using MFront} A user-defined material created with MFront can be used if MGIS is installed \newline diff --git a/src/sm/CMakeLists.txt b/src/sm/CMakeLists.txt index 9f2bd8a72..bfa78524b 100644 --- a/src/sm/CMakeLists.txt +++ b/src/sm/CMakeLists.txt @@ -97,6 +97,7 @@ set (sm_element Elements/LatticeElements/latticebeam3d.C Elements/LatticeElements/latticebeam3dboundary.C Elements/LatticeElements/latticeframe3d.C + Elements/LatticeElements/latticeframe3dnl.C Elements/tet21ghostsolid.C Elements/quad1platesubsoil.C Elements/quad2platesubsoil.C @@ -297,9 +298,9 @@ set (sm_material Materials/LatticeMaterials/latticeplasticitydamageviscoelastic.C Materials/LatticeMaterials/latticeframeelastic.C Materials/LatticeMaterials/latticeframesteelplastic.C + Materials/LatticeMaterials/latticeframeconcreteplastic.C Materials/1D/isoasymm1d.C ) - if (USE_MFRONT) list (APPEND sm_material Materials/mfrontusermaterial.C) endif () @@ -429,8 +430,9 @@ set (sm ${sm_obsolete} ${sm_new} ${sm_boundary_conditions} - ${sm_quasicontinuum} - ) + ${sm_quasicontinuum} + Elements/LatticeElements/latticebeam3dboundary.h +) if (USE_PARALLEL) list (APPEND sm ${sm_parallel}) diff --git a/src/sm/Elements/LatticeElements/latticeframe3d.C b/src/sm/Elements/LatticeElements/latticeframe3d.C index 796007989..542a7eaa0 100644 --- a/src/sm/Elements/LatticeElements/latticeframe3d.C +++ b/src/sm/Elements/LatticeElements/latticeframe3d.C @@ -101,14 +101,14 @@ LatticeFrame3d::computeBmatrixAt(GaussPoint *aGaussPoint, FloatMatrix &answer, i answer.at(2, 3) = 0.; answer.at(2, 4) = 0.; answer.at(2, 5) = 0; - answer.at(2, 6) = -this->length*(1-this->s)/2.; + answer.at(2, 6) = -this->length*(1.-this->s)/2.; //Second node answer.at(2, 7) = 0.; answer.at(2, 8) = 1.; answer.at(2, 9) = 0.; answer.at(2, 10) = 0.; answer.at(2, 11) = 0; - answer.at(2, 12) = -this->length*(1+this->s)/2.; + answer.at(2, 12) = -this->length*(1.+this->s)/2.; //Shear displacement jump in z-plane //first node @@ -116,14 +116,14 @@ LatticeFrame3d::computeBmatrixAt(GaussPoint *aGaussPoint, FloatMatrix &answer, i answer.at(3, 2) = 0.; answer.at(3, 3) = -1.; answer.at(3, 4) = 0.; - answer.at(3, 5) = this->length*(1-this->s)/2.; + answer.at(3, 5) = this->length*(1.-this->s)/2.; answer.at(3, 6) = 0.; //Second node answer.at(3, 7) = 0.; answer.at(3, 8) = 0.; answer.at(3, 9) = 1.; answer.at(3, 10) = 0.; - answer.at(3, 11) = this->length*(1+this->s)/2.; + answer.at(3, 11) = this->length*(1.+this->s)/2.; answer.at(3, 12) = 0.; //Rotation around x-axis @@ -241,7 +241,8 @@ LatticeFrame3d::computeStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMod dbj.times(1. / length); bjt.beTranspositionOf(bj); answer.beProductOf(bjt, dbj); - + //printf("answer/n"); + //answer.printYourself(); return; } @@ -305,6 +306,7 @@ LatticeFrame3d::giveInternalForcesVector(FloatArray &answer, answer.clear(); this->computeBmatrixAt(integrationRulesArray [ 0 ]->getIntegrationPoint(0), b); + bt.beTranspositionOf(b); if ( useUpdatedGpRecord == 1 ) { @@ -318,17 +320,13 @@ LatticeFrame3d::giveInternalForcesVector(FloatArray &answer, strain.times(1./this->length); this->computeStressVector(stress, strain, integrationRulesArray [ 0 ]->getIntegrationPoint(0), tStep); } - answer.beProductOf(bt, stress); - if ( !this->isActivated(tStep) ) { answer.zero(); return; } } - - bool LatticeFrame3d::computeGtoLRotationMatrix(FloatMatrix &answer) { @@ -357,10 +355,10 @@ LatticeFrame3d::giveLocalCoordinateSystem(FloatMatrix &answer) Node *nodeA, *nodeB; nodeA = this->giveNode(1); nodeB = this->giveNode(2); - lx.beDifferenceOf(nodeB->giveCoordinates(), nodeA->giveCoordinates() ); lx.normalize(); + if ( this->referenceNode ) { Node *refNode = this->giveDomain()->giveNode(this->referenceNode); help.beDifferenceOf(refNode->giveCoordinates(), nodeA->giveCoordinates() ); @@ -446,6 +444,7 @@ LatticeFrame3d::initializeFrom(InputRecord &ir) this->s = 0.; IR_GIVE_OPTIONAL_FIELD(ir, s, _IFT_LatticeFrame3d_s); + } @@ -467,8 +466,6 @@ LatticeFrame3d::computeLength() return length; } - - void LatticeFrame3d::computeLumpedMassMatrix(FloatMatrix &answer, TimeStep *tStep) // Returns the lumped mass matrix of the receiver. This expression is diff --git a/src/sm/Elements/LatticeElements/latticeframe3d.h b/src/sm/Elements/LatticeElements/latticeframe3d.h index e7e9fbf81..c22c8dda0 100644 --- a/src/sm/Elements/LatticeElements/latticeframe3d.h +++ b/src/sm/Elements/LatticeElements/latticeframe3d.h @@ -110,7 +110,7 @@ class LatticeFrame3d : public LatticeStructuralElement protected: - void computeBmatrixAt(GaussPoint *, FloatMatrix &, int = 1, int = ALL_STRAINS) override; + virtual void computeBmatrixAt(GaussPoint *, FloatMatrix &, int = 1, int = ALL_STRAINS) override; bool computeGtoLRotationMatrix(FloatMatrix &) override; void computeLumpedMassMatrix(FloatMatrix &answer, TimeStep *tStep) override; void computeMassMatrix(FloatMatrix &answer, TimeStep *tStep) override diff --git a/src/sm/Elements/LatticeElements/latticeframe3dnl.C b/src/sm/Elements/LatticeElements/latticeframe3dnl.C new file mode 100755 index 000000000..63a5bf950 --- /dev/null +++ b/src/sm/Elements/LatticeElements/latticeframe3dnl.C @@ -0,0 +1,619 @@ +/* + * + * ##### ##### ###### ###### ### ### + * ## ## ## ## ## ## ## ### ## + * ## ## ## ## #### #### ## # ## + * ## ## ## ## ## ## ## ## + * ## ## ## ## ## ## ## ## + * ##### ##### ## ###### ## ## + * + * + * OOFEM : Object Oriented Finite Element Code + * + * Copyright (C) 1993 - 2023 Borek Patzak + * + * + * + * Czech Technical University, Faculty of Civil Engineering, + * Department of Structural Mechanics, 166 29 Prague, Czech Republic + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + */ + +#include "domain.h" +#include "latticeframe3dnl.h" +#include "../sm/Materials/LatticeMaterials/latticematstatus.h" +#include "node.h" +#include "material.h" +#include "gausspoint.h" +#include "gaussintegrationrule.h" +#include "floatmatrix.h" +#include "floatmatrixf.h" +#include "intarray.h" +#include "floatarray.h" +#include "floatarrayf.h" +#include "mathfem.h" +#include "latticeframe3d.h" +#include "contextioerr.h" +#include "datastream.h" +#include "classfactory.h" +#include "sm/CrossSections/latticecrosssection.h" +#include "engngm.h" + +#ifdef __OOFEG +#include "oofeggraphiccontext.h" + #include "../sm/Materials/structuralmaterial.h" +#endif + +namespace oofem { + REGISTER_Element(LatticeFrame3dNL); + + LatticeFrame3dNL::LatticeFrame3dNL(int n, Domain *aDomain) : LatticeFrame3d(n, aDomain) + { + numberOfDofMans = 2; + } + + LatticeFrame3dNL::~LatticeFrame3dNL() + {} + void + LatticeFrame3dNL::computeBmatrixAt(GaussPoint *aGaussPoint, FloatMatrix &answer, int li, int ui) + //Returns the strain matrix of the receiver. + { + FloatArray u; + TimeStep *tStep = this->domain->giveEngngModel()->giveCurrentStep(); + this->computeVectorOf(VM_Total, tStep, u); + + this->length = computeLength(); + double l1 = this->length * ( 1. - this->s ) / 2; + double l2 = this->length * ( 1. + this->s ) / 2; + answer.resize(6, 12); + answer.zero(); +// double cx1=(cos(u.at(5))*cos(u.at(6)))*l1; +// double cy1=(cos(u.at(4))*sin(u.at(6))+sin(u.at(4))*sin(u.at(5))*cos(u.at(6)))*l1; +// double cz1=(sin(u.at(4))*sin(u.at(6))-cos(u.at(4))*sin(u.at(5))*cos(u.at(6)))*l1; +// +// +// double cx2=(cos(u.at(11))*cos(u.at(12)))*l2; +// double cy2=(cos(u.at(10))*sin(u.at(12))+sin(u.at(10))*sin(u.at(11))*cos(u.at(12)))*l2; +// double cz2=(sin(u.at(10))*sin(u.at(12))-cos(u.at(10))*sin(u.at(11))*cos(u.at(12)))*l2; + + double cx1=(cos(u.at(5))*cos(u.at(6)))*l1; + double cy1=(cos(u.at(5))*sin(u.at(6)))*l1; + double cz1=(-sin(u.at(5)))*l1; + + + double cx2=(cos(u.at(11))*cos(u.at(12)))*l2; + double cy2=(cos(u.at(11))*sin(u.at(12)))*l2; + double cz2=(-sin(u.at(11)))*l2; + + + //Normal displacement jump in x-direction + //First node + answer.at(1, 1) = -1.; + answer.at(1, 2) = 0.; + answer.at(1, 3) = 0.; + answer.at(1, 4) = 0.; + answer.at(1, 5) = -cz1; + answer.at(1, 6) = cy1; + //Second node + answer.at(1, 7) = 1.; + answer.at(1, 8) = 0.; + answer.at(1, 9) = 0.; + answer.at(1, 10) = 0.; + answer.at(1, 11) = -cz2; + answer.at(1, 12) = cy2; + + //Shear displacement jump in y-plane + //first node + answer.at(2, 1) = 0.; + answer.at(2, 2) = -1.; + answer.at(2, 3) = 0.; + answer.at(2, 4) = cz1; + answer.at(2, 5) = 0; + answer.at(2, 6) = -cx1; + //Second node + answer.at(2, 7) = 0.; + answer.at(2, 8) = 1.; + answer.at(2, 9) = 0.; + answer.at(2, 10) = cz2; + answer.at(2, 11) = 0; + answer.at(2, 12) = -cx2; + + //Shear displacement jump in z-plane + //first node + answer.at(3, 1) = 0.; + answer.at(3, 2) = 0.; + answer.at(3, 3) = -1.; + answer.at(3, 4) = -cy1; + answer.at(3, 5) = cx1; + answer.at(3, 6) = 0.; + //Second node + answer.at(3, 7) = 0.; + answer.at(3, 8) = 0.; + answer.at(3, 9) = 1.; + answer.at(3, 10) = -cy2; + answer.at(3, 11) = cx2; + answer.at(3, 12) = 0.; + + //Rotation around x-axis + //First node + answer.at(4, 1) = 0.; + answer.at(4, 2) = 0; + answer.at(4, 3) = 0.; + answer.at(4, 4) = -1.; + answer.at(4, 5) = 0.; + answer.at(4, 6) = 0.; + //Second node + answer.at(4, 7) = 0.; + answer.at(4, 8) = 0.; + answer.at(4, 9) = 0.; + answer.at(4, 10) = 1.; + answer.at(4, 11) = 0.; + answer.at(4, 12) = 0.; + + //Rotation around y-axis + //First node + answer.at(5, 1) = 0.; + answer.at(5, 2) = 0.; + answer.at(5, 3) = 0.; + answer.at(5, 4) = 0.; + answer.at(5, 5) = -1.; + answer.at(5, 6) = 0.; + //Second node + answer.at(5, 7) = 0.; + answer.at(5, 8) = 0.; + answer.at(5, 9) = 0.; + answer.at(5, 10) = 0.; + answer.at(5, 11) = 1.; + answer.at(5, 12) = 0.; + + //Rotation around z-axis + //First node + answer.at(6, 1) = 0.; + answer.at(6, 2) = 0.; + answer.at(6, 3) = 0.; + answer.at(6, 4) = 0.; + answer.at(6, 5) = 0.; + answer.at(6, 6) = -1.; + //Second node + answer.at(6, 7) = 0.; + answer.at(6, 8) = 0.; + answer.at(6, 9) = 0.; + answer.at(6, 10) = 0.; + answer.at(6, 11) = 0.; + answer.at(6, 12) = 1.; +// + return; + } + void + LatticeFrame3dNL::computeStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, + TimeStep *tStep) + { + FloatMatrix d, bt, db, b; + FloatArray u; + + this->computeVectorOf(VM_Total, tStep, u); + this->length = computeLength(); + + answer.resize(12, 12); + answer.zero(); + this->computeBmatrixAt(integrationRulesArray [ 0 ]->getIntegrationPoint(0), b); + this->computeConstitutiveMatrixAt(d, rMode, integrationRulesArray [ 0 ]->getIntegrationPoint(0), tStep); + + db.beProductOf(d, b); + db.times(1. / length); + bt.beTranspositionOf(b); + answer.beProductOf(bt, db); + +// printf("answer/n"); +// answer.printYourself(); + return; + } +// void +// LatticeFrame3dNL::computeStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, +// TimeStep *tStep) +// { +// FloatMatrix d; +// FloatArray u; +// this->length = computeLength(); +// answer.resize(12, 12); +// answer.zero(); +// +// this->computeConstitutiveMatrixAt(d, rMode, integrationRulesArray [ 0 ]->getIntegrationPoint(0), tStep); +// +// double l1 = this->length * ( 1. - this->s ) / 2; +// double l2 = this->length * ( 1. + this->s ) / 2; +// +// this->computeVectorOf(VM_Total, tStep, u); +// +// double cx1=(cos(u.at(5))*cos(u.at(6)))*l1; +// double cy1=(cos(u.at(4))*sin(u.at(6))+sin(u.at(4))*sin(u.at(5))*cos(u.at(6)))*l1; +// double cz1=(sin(u.at(4))*sin(u.at(6))-cos(u.at(4))*sin(u.at(5))*cos(u.at(6)))*l1; +// +// +// double cx2=(cos(u.at(11))*cos(u.at(12)))*l2; +// double cy2=(cos(u.at(10))*sin(u.at(12))+sin(u.at(10))*sin(u.at(11))*cos(u.at(12)))*l2; +// double cz2=(sin(u.at(10))*sin(u.at(12))-cos(u.at(10))*sin(u.at(11))*cos(u.at(12)))*l2; +// +// +// +// //Axial 1 +// answer.at(1, 1) = d.at(1, 1); +// answer.at(1, 2) = 0.; +// answer.at(1, 3) = 0.; +// answer.at(1, 4) = 0.; +// answer.at(1, 5) = -d.at(1,1)*(-cz1); +// answer.at(1, 6) = -d.at(1,1)*(cy1); +// answer.at(1, 7) = -d.at(1, 1); +// answer.at(1, 8) = 0.; +// answer.at(1, 9) = 0.; +// answer.at(1, 10) = 0.; +// answer.at(1, 11) = -d.at(1,1)*(-cz2); +// answer.at(1, 12) = -d.at(1,1)*(cy2); +// +// //Shear Y 1 +// answer.at(2, 1) = 0; +// answer.at(2, 2) = d.at(2, 2); +// answer.at(2, 3) = 0.; +// answer.at(2, 4) = -d.at(2,2)*(cz1); +// answer.at(2, 5) = 0; +// answer.at(2, 6) = -d.at(2,2)*(-cx1); +// answer.at(2, 7) = 0.; +// answer.at(2, 8) = -d.at(2,2); +// answer.at(2, 9) = 0.; +// answer.at(2, 10) = -d.at(2,2)*(cz2); +// answer.at(2, 11) = 0; +// answer.at(2, 12) = -d.at(2,2)*(-cx2); +// +// //Shear Z 1 +// answer.at(3, 1) = 0; +// answer.at(3, 2) = 0.; +// answer.at(3, 3) = d.at(3, 3); +// answer.at(3, 4) = -d.at(3,3)*(-cy1); +// answer.at(3, 5) = -d.at(3,3)*(cx1); +// answer.at(3, 6) = 0; +// answer.at(3, 7) = 0.; +// answer.at(3, 8) = 0; +// answer.at(3, 9) = -d.at(3, 3); +// answer.at(3, 10) = -d.at(3,3)*(-cy2); +// answer.at(3, 11) = -d.at(3,3)*(cx2); +// answer.at(3, 12) = 0; +// +// // Mx 1 +// answer.at(4, 1) = 0; +// answer.at(4, 2) = d.at(2,2)*(-cz1); +// answer.at(4, 3) = -d.at(3,3)*(-cy1); +// answer.at(4, 4) = -d.at(3,3)*(-cy1*cy1)+d.at(2,2)*(cz1*cz1)+d.at(4,4); +// answer.at(4, 5) = -d.at(3,3)*(cy1*cx1); +// answer.at(4, 6) = d.at(2,2)*(-cz1*cx1); +// answer.at(4, 7) = 0.; +// answer.at(4, 8) = d.at(2,2)*(cz1); +// answer.at(4, 9) = -d.at(3,3)*(cy1); +// answer.at(4, 10) = d.at(3,3)*(cy1*cy2)+d.at(2,2)*(+cz1*cz2)-d.at(4,4); +// answer.at(4, 11) = -d.at(3,3)*(cy1*cx2); +// answer.at(4, 12) = d.at(2,2)*(-cz1*cx2); +// +// // My 1 +// answer.at(5, 1) = -d.at(1,1)*(-cz1); +// answer.at(5, 2) = 0; +// answer.at(5, 3) = d.at(3,3)*(-cx1); +// answer.at(5, 4) = d.at(3,3)*(-cx1*cy1); +// answer.at(5, 5) = -d.at(1,1)*(-cz1*cz1)+d.at(3,3)*(cx1*cx1)+d.at(5,5); +// answer.at(5, 6) = -d.at(1,1)*(+cz1*cy1); +// answer.at(5, 7) = -d.at(1,1)*(cz1); +// answer.at(5, 8) = 0; +// answer.at(5, 9) = d.at(3,3)*(cx1); +// answer.at(5, 10) = d.at(3,3)*(-cx1*cy2); +// answer.at(5, 11) = -d.at(1,1)*(-cz1*cz2)+d.at(3,3)*(cx1*cx2)-d.at(5,5); +// answer.at(5, 12) = -d.at(1,1)*(cz1*cy2); +// +// // Mz 1 +// answer.at(6, 1) = d.at(1,1)*(-cy1); +// answer.at(6, 2) = -d.at(2,2)*(-cx1); +// answer.at(6, 3) = 0; +// answer.at(6, 4) = -d.at(2,2)*(cx1*cz1); +// answer.at(6, 5) = d.at(1,1)*(-cy1*cz1); +// answer.at(6, 6) = d.at(1,1)*(cy1*cy1)-d.at(2,2)*(-cx1*cx1)+d.at(6,6); +// answer.at(6, 7) = d.at(1,1)*(cy1); +// answer.at(6, 8) = -d.at(2,2)*(cx1); +// answer.at(6, 9) = 0; +// answer.at(6, 10) = -d.at(2,2)*(cx1*cz2); +// answer.at(6, 11) = d.at(1,1)*(-cy1*cz2); +// answer.at(6, 12) = d.at(1,1)*(cy1*cy2)-d.at(2,2)*(-cx1*cx2)-d.at(6,6); +// +// //Axial 2 +// answer.at(7, 1) = -d.at(1, 1); +// answer.at(7, 2) = 0.; +// answer.at(7, 3) = 0.; +// answer.at(7, 4) = 0.; +// answer.at(7, 5) = d.at(1,1)*(-cz1); +// answer.at(7, 6) = d.at(1,1)*(cy1); +// answer.at(7, 7) = d.at(1, 1); +// answer.at(7, 8) = 0.; +// answer.at(7, 9) = 0.; +// answer.at(7, 10) = 0.; +// answer.at(7, 11) = d.at(1,1)*(-cz2); +// answer.at(7, 12) = d.at(1,1)*(cy2); +// +// //Shear Y 2 +// answer.at(8, 1) = 0; +// answer.at(8, 2) = -d.at(2, 2); +// answer.at(8, 3) = 0.; +// answer.at(8, 4) = d.at(2,2)*(cz1); +// answer.at(8, 5) = 0; +// answer.at(8, 6) = d.at(2,2)*(-cx1); +// answer.at(8, 7) = 0.; +// answer.at(8, 8) = d.at(2,2); +// answer.at(8, 9) = 0.; +// answer.at(8, 10) = d.at(2,2)*(cz2); +// answer.at(8, 11) = 0; +// answer.at(8, 12) = d.at(2,2)*(-cx2); +// +// //Shear Z 2 +// answer.at(9, 1) = 0; +// answer.at(9, 2) = 0.; +// answer.at(9, 3) = -d.at(3, 3); +// answer.at(9, 4) = d.at(3,3)*(-cy1); +// answer.at(9, 5) = d.at(3,3)*(cx1); +// answer.at(9, 6) = 0; +// answer.at(9, 7) = 0.; +// answer.at(9, 8) = 0; +// answer.at(9, 9) = d.at(3, 3); +// answer.at(9, 10) = d.at(3,3)*(-cy2); +// answer.at(9, 11) = d.at(3,3)*(cx2); +// answer.at(9, 12) = 0; +// +// // Mx 2 +// answer.at(10, 1) = 0; +// answer.at(10, 2) = -d.at(2,2)*(cz2); +// answer.at(10, 3) = d.at(3,3)*(cy2); +// answer.at(10, 4) = d.at(3,3)*(cy2*cy1)-d.at(2,2)*(-cz2*cz1)-d.at(4,4); +// answer.at(10, 5) = d.at(3,3)*(-cy2*cx1); +// answer.at(10, 6) = -d.at(2,2)*(cz2*cx1); +// answer.at(10, 7) = 0.; +// answer.at(10, 8) = -d.at(2,2)*(-cz2); +// answer.at(10, 9) = d.at(3,3)*(-cy2); +// answer.at(10, 10) = d.at(3,3)*(cy2*cy2)-d.at(2,2)*(-cz2*cz2)+d.at(4,4); +// answer.at(10, 11) = d.at(3,3)*(-cy2*cx2); +// answer.at(10, 12) = -d.at(2,2)*(+cz2*cx2); +// +// // My 2 +// answer.at(11, 1) = d.at(1,1)*(cz2); +// answer.at(11, 2) = 0; +// answer.at(11, 3) = -d.at(3,3)*(cx2); +// answer.at(11, 4) = -d.at(3,3)*(cx2*cy1); +// answer.at(11, 5) = d.at(1,1)*(cz2*cz1)-d.at(3,3)*(-cx2*cx1)-d.at(5,5); +// answer.at(11, 6) = d.at(1,1)*(-cz2*cy1); +// answer.at(11, 7) = d.at(1,1)*(-cz2); +// answer.at(11, 8) = 0; +// answer.at(11, 9) = -d.at(3,3)*(-cx2); +// answer.at(11, 10) = -d.at(3,3)*(cx2*cy2); +// answer.at(11, 11) = d.at(1,1)*(cz2*cz2)-d.at(3,3)*(-cx2*cx2)+d.at(5,5); +// answer.at(11, 12) = d.at(1,1)*(-cz2*cy2); +// +// // Mz 2 +// answer.at(12, 1) = -cy2*d.at(1,1); +// answer.at(12, 2) = cx2*d.at(2,2); +// answer.at(12, 3) = 0; +// answer.at(12, 4) = -cx2*d.at(2,2)*(cz1); +// answer.at(12, 5) = cy2*d.at(1,1)*(-cz1); +// answer.at(12, 6) = cy2*d.at(1,1)*(cy1)-cx2*d.at(2,2)*(-cx1)-d.at(6,6); +// answer.at(12, 7) = cy2*d.at(1,1); +// answer.at(12, 8) = -cx2*d.at(2,2); +// answer.at(12, 9) = 0; +// answer.at(12, 10) = -cx2*d.at(2,2)*(cz2); +// answer.at(12, 11) = cy2*d.at(1,1)*(-cz2); +// answer.at(12, 12) = cy2*d.at(1,1)*(cy2)-cx2*d.at(2,2)*(-cx2)+d.at(6,6); +// +// answer.times(1. / this->length); +// return; +// } + + void + LatticeFrame3dNL::computeStrainVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep) + // Computes the vector containing the strains at the Gauss point gp of + // the receiver, at time step tStep. The nature of these strains depends + // on the element's type. + { + FloatArray u; + this->computeVectorOf(VM_Total, tStep, u); + +// printf("displacement n:"); +// un.printYourself(); + + this->length = computeLength(); + double l1 = this->length * ( 1. - this->s ) / 2; + double l2 = this->length * ( 1. + this->s ) / 2; + // double ln = sqrt(pow(njn.at(1)-nin.at(1), 2) + pow(njn.at(2)-nin.at(2), 2) + pow(njn.at(3)-nin.at(3), 2) ); + LatticeMaterialStatus *lmatStat = dynamic_cast < LatticeMaterialStatus * > ( integrationRulesArray [ 0 ]->getIntegrationPoint(0)->giveMaterialStatus() ); + auto strain = lmatStat->giveLatticeStrain(); + +// double cx1=(cos(u.at(5))*cos(u.at(6)))*l1; +// double cy1=(cos(u.at(4))*sin(u.at(6))+sin(u.at(4))*sin(u.at(5))*cos(u.at(6)))*l1; +// double cz1=(sin(u.at(4))*sin(u.at(6))-cos(u.at(4))*sin(u.at(5))*cos(u.at(6)))*l1; +// +// +// double cx2=(cos(u.at(11))*cos(u.at(12)))*l2; +// double cy2=(cos(u.at(10))*sin(u.at(12))+sin(u.at(10))*sin(u.at(11))*cos(u.at(12)))*l2; +// double cz2=(sin(u.at(10))*sin(u.at(12))-cos(u.at(10))*sin(u.at(11))*cos(u.at(12)))*l2; + double cx1=(cos(u.at(5))*cos(u.at(6)))*l1; + double cy1=(cos(u.at(5))*sin(u.at(6)))*l1; + double cz1=(-sin(u.at(5)))*l1; + + + double cx2=(cos(u.at(11))*cos(u.at(12)))*l2; + double cy2=(cos(u.at(11))*sin(u.at(12)))*l2; + double cz2=(-sin(u.at(11)))*l2; + // + answer.resize(6); + answer.at(1) = u.at(7)-u.at(1)-cx2 -cx1 +l1+l2; + answer.at(2) = u.at(8)-u.at(2)-cy2 -cy1; + answer.at(3) = u.at(9)-u.at(3)-cz2 -cz1; + answer.at(4) = u.at(10)-u.at(4); + answer.at(5) = u.at(11)-u.at(5); + answer.at(6) = u.at(12)-u.at(6); + answer.times(1. / this->length); + + } + // + void + LatticeFrame3dNL::giveInternalForcesVector(FloatArray &answer, + TimeStep *tStep, int useUpdatedGpRecord) + { + FloatMatrix b, bt, bf, d; + FloatArray u, stress, strain; + this->computeVectorOf(VM_Total, tStep, u); + // this->computeVectorOf(VM_Incremental, tStep, u); + this->length = computeLength(); + GaussPoint *gp = this->integrationRulesArray [ 0 ]->getIntegrationPoint(0); + // Total stress + this->computeStrainVector(strain, gp, tStep); + this->computeStressVector(stress, strain, integrationRulesArray [ 0 ]->getIntegrationPoint(0), tStep); + // Old stresses + LatticeMaterialStatus *lmatStat = dynamic_cast < LatticeMaterialStatus * > ( integrationRulesArray [ 0 ]->getIntegrationPoint(0)->giveMaterialStatus() ); + auto oldStress = lmatStat->giveLatticeStress(); + // + auto oldInternalForces = lmatStat->giveInternalForces(); + double l1 = this->length * ( 1. - this->s ) / 2; + double l2 = this->length * ( 1. + this->s ) / 2; + // +// double cx1=(cos(u.at(5))*cos(u.at(6)))*l1; +// double cy1=(cos(u.at(4))*sin(u.at(6))+sin(u.at(4))*sin(u.at(5))*cos(u.at(6)))*l1; +// double cz1=(sin(u.at(4))*sin(u.at(6))-cos(u.at(4))*sin(u.at(5))*cos(u.at(6)))*l1; +// +// +// double cx2=(cos(u.at(11))*cos(u.at(12)))*l2; +// double cy2=(cos(u.at(10))*sin(u.at(12))+sin(u.at(10))*sin(u.at(11))*cos(u.at(12)))*l2; +// double cz2=(sin(u.at(10))*sin(u.at(12))-cos(u.at(10))*sin(u.at(11))*cos(u.at(12)))*l2; + + double cx1=(cos(u.at(5))*cos(u.at(6)))*l1; + double cy1=(cos(u.at(5))*sin(u.at(6)))*l1; + double cz1=(-sin(u.at(5)))*l1; + + + double cx2=(cos(u.at(11))*cos(u.at(12)))*l2; + double cy2=(cos(u.at(11))*sin(u.at(12)))*l2; + double cz2=(-sin(u.at(11)))*l2; + // + // + answer.resize(12); + answer.at(1) = -stress.at(1); + answer.at(2) = -stress.at(2); + answer.at(3) = -stress.at(3); + answer.at(4) = stress.at(2)*cz1- stress.at(3)*cy1- stress.at(4); + answer.at(5) = -stress.at(1)*cz1+ stress.at(3)*cx1- stress.at(5); + answer.at(6) = stress.at(1)*cy1- stress.at(2)*cx1- stress.at(6); + answer.at(7) = stress.at(1); + answer.at(8) = stress.at(2); + answer.at(9) = stress.at(3); + answer.at(10) = stress.at(2)*cz2- stress.at(3)*cy2+ stress.at(4); + answer.at(11) = -stress.at(1)*cz2+ stress.at(3)*cx2+ stress.at(5); + answer.at(12) = stress.at(1)*cy2- stress.at(2)*cx2+ stress.at(6); + // + lmatStat->letTempInternalForcesBe(answer); + + } + + int + LatticeFrame3dNL::giveLocalCoordinateSystem(FloatMatrix &answer) + { + FloatArray lx, ly, lz, help(3); + FloatArray coordA, coordB; + FloatArray uA(6), uAIncr(6), uB(6), uBIncr(6); + IntArray dofid = { + 1, 2, 3, 4, 5, 6 + }; + + TimeStep *tStep = this->domain->giveEngngModel()->giveCurrentStep(); + + Node *nodeA, *nodeB; + nodeA = this->giveNode(1); + nodeB = this->giveNode(2); + + //Local coordinate system is determined from the displacement of last step. + + coordA = nodeA->giveCoordinates(); + nodeA->giveUnknownVector(uA, dofid, VM_Total, tStep, false); + nodeA->giveUnknownVector(uAIncr, dofid, VM_Incremental, tStep, false); +// for (int i = 1; i <= 3; i++) { +// coordA.at(i) += uA.at(i) - uAIncr.at(i); +// } + + coordB = nodeB->giveCoordinates(); + nodeB->giveUnknownVector(uB, dofid, VM_Total, tStep, false); + nodeB->giveUnknownVector(uBIncr, dofid, VM_Incremental, tStep, false); + +// for (int i = 1; i <= 3; i++) { +// coordB.at(i) += uB.at(i) - uBIncr.at(i); +// } + + lx.beDifferenceOf(coordB, coordA); + lx.normalize(); + + + if ( this->referenceNode ) { + Node *refNode = this->giveDomain()->giveNode(this->referenceNode); + help.beDifferenceOf( refNode->giveCoordinates(), nodeA->giveCoordinates() ); + + lz.beVectorProductOf(lx, help); + lz.normalize(); + } else if ( this->zaxis.giveSize() > 0 ) { + lz = this->zaxis; + lz.add(lz.dotProduct(lx), lx); + lz.normalize(); + } else { + FloatMatrix rot(3, 3); + double theta = referenceAngle * M_PI / 180.0; + + rot.at(1, 1) = cos(theta) + pow(lx.at(1), 2) * ( 1 - cos(theta) ); + rot.at(1, 2) = lx.at(1) * lx.at(2) * ( 1 - cos(theta) ) - lx.at(3) * sin(theta); + rot.at(1, 3) = lx.at(1) * lx.at(3) * ( 1 - cos(theta) ) + lx.at(2) * sin(theta); + + rot.at(2, 1) = lx.at(2) * lx.at(1) * ( 1 - cos(theta) ) + lx.at(3) * sin(theta); + rot.at(2, 2) = cos(theta) + pow(lx.at(2), 2) * ( 1 - cos(theta) ); + rot.at(2, 3) = lx.at(2) * lx.at(3) * ( 1 - cos(theta) ) - lx.at(1) * sin(theta); + + rot.at(3, 1) = lx.at(3) * lx.at(1) * ( 1 - cos(theta) ) - lx.at(2) * sin(theta); + rot.at(3, 2) = lx.at(3) * lx.at(2) * ( 1 - cos(theta) ) + lx.at(1) * sin(theta); + rot.at(3, 3) = cos(theta) + pow(lx.at(3), 2) * ( 1 - cos(theta) ); + + help.at(3) = 1.0; // up-vector + + // here is ly is used as a temp var + if ( fabs( lx.dotProduct(help) ) > 0.999 ) { // Check if it is vertical + ly = { + 0., 1., 0. + }; + } else { + ly.beVectorProductOf(lx, help); + } + lz.beProductOf(rot, ly); + lz.normalize(); + } + + ly.beVectorProductOf(lz, lx); + ly.normalize(); + + answer.resize(3, 3); + answer.zero(); + for ( int i = 1; i <= 3; i++ ) { + answer.at(1, i) = lx.at(i); + answer.at(2, i) = ly.at(i); + answer.at(3, i) = lz.at(i); + } + + return 1; + } +} // end namespace oofem \ No newline at end of file diff --git a/src/sm/Elements/LatticeElements/latticeframe3dnl.h b/src/sm/Elements/LatticeElements/latticeframe3dnl.h new file mode 100644 index 000000000..dcc9e692f --- /dev/null +++ b/src/sm/Elements/LatticeElements/latticeframe3dnl.h @@ -0,0 +1,76 @@ +/* + * + * ##### ##### ###### ###### ### ### + * ## ## ## ## ## ## ## ### ## + * ## ## ## ## #### #### ## # ## + * ## ## ## ## ## ## ## ## + * ## ## ## ## ## ## ## ## + * ##### ##### ## ###### ## ## + * + * + * OOFEM : Object Oriented Finite Element Code + * + * Copyright (C) 1993 - 2023 Borek Patzak + * + * + * + * Czech Technical University, Faculty of Civil Engineering, + * Department of Structural Mechanics, 166 29 Prague, Czech Republic + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + */ + +#ifndef latticeframe3dnl_h +#define latticeframe3dnl_h + +#include "latticeframe3d.h" + +///@name Input fields for LatticeFrame3dNL +//@{ +#define _IFT_LatticeFrame3dNL_Name "latticeframe3dnl" + +//@} + +namespace oofem { +/** +This class implements a geometric nonlinear 3-dimensional frame element. It is an extension of a geometric linear 3-dimensional frame element based on rigid body spring theory called LatticeFrame3D presented in Toi 1991 and Toi 1993. It belongs to the group of lattice models in the OOFEM structure, but can be used as a standard 3D beam element. +References: +Toi, Y. (1991). Shifted integration technique in one‐dimensional plastic collapse analysis using linear and cubic finite elements. International Journal for Numerical Methods in Engineering, 31(8), 1537-1552. +Toi, Y., & Isobe, D. (1993). Adaptively shifted integration technique for finite element collapse analysis of framed structures. International Journal for Numerical Methods in Engineering, 36(14), 2323-2339. +Authors: Gumaa Abdelrhim and Peter Grassl, 2023 +*/ + + class LatticeFrame3dNL : public LatticeFrame3d + { + protected: + + public: + LatticeFrame3dNL(int n, Domain *); + virtual ~LatticeFrame3dNL(); + + int giveLocalCoordinateSystem(FloatMatrix &answer) override; + + const char *giveInputRecordName() const override { return _IFT_LatticeFrame3dNL_Name; } + const char *giveClassName() const override { return "latticeframe3dnl"; } + + void giveInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord) override; + + protected: + virtual void computeBmatrixAt(GaussPoint *, FloatMatrix &, int = 1, int = ALL_STRAINS); + void computeStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep) override; + virtual void computeStrainVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep) override; + }; +} // end namespace oofem +#endif \ No newline at end of file diff --git a/src/sm/Elements/LatticeElements/latticeframe3dnl2.C b/src/sm/Elements/LatticeElements/latticeframe3dnl2.C new file mode 100644 index 000000000..be93ea846 --- /dev/null +++ b/src/sm/Elements/LatticeElements/latticeframe3dnl2.C @@ -0,0 +1,493 @@ +/* + * + * ##### ##### ###### ###### ### ### + * ## ## ## ## ## ## ## ### ## + * ## ## ## ## #### #### ## # ## + * ## ## ## ## ## ## ## ## + * ## ## ## ## ## ## ## ## + * ##### ##### ## ###### ## ## + * + * + * OOFEM : Object Oriented Finite Element Code + * + * Copyright (C) 1993 - 2023 Borek Patzak + * + * + * + * Czech Technical University, Faculty of Civil Engineering, + * Department of Structural Mechanics, 166 29 Prague, Czech Republic + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + */ + +#include "domain.h" +#include "latticeframe3dnl2.h" +#include "../sm/Materials/LatticeMaterials/latticematstatus.h" +#include "node.h" +#include "material.h" +#include "gausspoint.h" +#include "gaussintegrationrule.h" +#include "floatmatrix.h" +#include "floatmatrixf.h" +#include "intarray.h" +#include "floatarray.h" +#include "floatarrayf.h" +#include "mathfem.h" +#include "latticeframe3d.h" +#include "contextioerr.h" +#include "datastream.h" +#include "classfactory.h" +#include "sm/CrossSections/latticecrosssection.h" +#include "engngm.h" + +#ifdef __OOFEG + #include "oofeggraphiccontext.h" + #include "../sm/Materials/structuralmaterial.h" +#endif + +namespace oofem { + REGISTER_Element(LatticeFrame3dNL2); + + LatticeFrame3dNL2::LatticeFrame3dNL2(int n, Domain *aDomain) : LatticeFrame3d(n, aDomain) + { + numberOfDofMans = 2; + } + + LatticeFrame3dNL2::~LatticeFrame3dNL2() + {} + void + LatticeFrame3dNL2::computeBmatrixAt(GaussPoint *aGaussPoint, FloatMatrix &answer, int li, int ui) + //Returns the strain matrix of the receiver. + { + FloatArray u; + TimeStep *tStep = this->domain->giveEngngModel()->giveCurrentStep(); + this->computeVectorOf(VM_Total, tStep, u); + + FloatArray coordA(3), coordB(3); + coordA = giveNode(1)->giveCoordinates(); + coordB = giveNode(2)->giveCoordinates(); + + this->length = computeLength(); + double l1 = this->length * ( 1. + this->s ) / 2.; + double l2 = this->length * ( 1. - this->s ) / 2.; + + double lx1 = ( coordB.at(1) - coordA.at(1) ) * ( 1. + this->s ) / 2.; + double ly1 = ( coordB.at(2) - coordA.at(2) ) * ( 1. + this->s ) / 2.; + double lz1 = ( coordB.at(3) - coordA.at(3) ) * ( 1. + this->s ) / 2.; + + double lx2 = ( coordB.at(1) - coordA.at(1) ) * ( 1. - this->s ) / 2.; + double ly2 = ( coordB.at(2) - coordA.at(2) ) * ( 1. - this->s ) / 2.; + double lz2 = ( coordB.at(3) - coordA.at(3) ) * ( 1. - this->s ) / 2.; + + answer.resize(6, 12); + answer.zero(); + + + double cx1Old = ( cos( u.at(5) ) * cos( u.at(6) ) ) * l1; + double cy1Old = ( cos( u.at(4) ) * sin( u.at(6) ) + sin( u.at(4) ) * sin( u.at(5) ) * cos( u.at(6) ) ) * l1; + double cz1Old = ( sin( u.at(4) ) * sin( u.at(6) ) - cos( u.at(4) ) * sin( u.at(5) ) * cos( u.at(6) ) ) * l1; + + + double cx2Old = ( cos( u.at(11) ) * cos( u.at(12) ) ) * l2; + double cy2Old = ( cos( u.at(10) ) * sin( u.at(12) ) + sin( u.at(10) ) * sin( u.at(11) ) * cos( u.at(12) ) ) * l2; + double cz2Old = ( sin( u.at(10) ) * sin( u.at(12) ) - cos( u.at(10) ) * sin( u.at(11) ) * cos( u.at(12) ) ) * l2; + + + //First node rotations + double thetaX1 = u.at(4); + double thetaY1 = u.at(5); + double thetaZ1 = u.at(6); + + //Second node rotations + double thetaX2 = u.at(10); + double thetaY2 = u.at(11); + double thetaZ2 = u.at(12); + + +// double cx1 = lz1 * sin(thetaY1) + lx1 * cos(thetaY1) * cos(thetaZ1) - ly1 * cos(thetaY1) * sin(thetaZ1); +// double cx2 = lz2 * sin(thetaY2) + lx2 * cos(thetaY2) * cos(thetaZ2) - ly2 * cos(thetaY2) * sin(thetaZ2); +// double cy1 = lx1 * ( cos(thetaX1) * sin(thetaZ1) + cos(thetaZ1) * sin(thetaX1) * sin(thetaY1) ) + ly1 * ( cos(thetaX1) * cos(thetaZ1) - sin(thetaX1) * sin(thetaY1) * sin(thetaZ1) ) - lz1 * cos(thetaY1) * sin(thetaX1); +// double cy2 = lx2 * ( cos(thetaX2) * sin(thetaZ2) + cos(thetaZ2) * sin(thetaX2) * sin(thetaY2) ) + ly2 * ( cos(thetaX2) * cos(thetaZ2) - sin(thetaX2) * sin(thetaY2) * sin(thetaZ2) ) - lz2 * cos(thetaY2) * sin(thetaX2); +// double cz1 = lx1 * ( sin(thetaX1) * sin(thetaZ1) - cos(thetaX1) * cos(thetaZ1) * sin(thetaY1) ) + ly1 * ( cos(thetaZ1) * sin(thetaX1) + cos(thetaX1) * sin(thetaY1) * sin(thetaZ1) ) + lz1 * cos(thetaX1) * cos(thetaY1); +// double cz2 = lx2 * ( sin(thetaX2) * sin(thetaZ2) - cos(thetaX2) * cos(thetaZ2) * sin(thetaY2) ) + ly2 * ( cos(thetaZ2) * sin(thetaX2) + cos(thetaX2) * sin(thetaY2) * sin(thetaZ2) ) + lz2 * cos(thetaX2) * cos(thetaY2); + + double cx1 = (cos(thetaZ1) * cos(thetaY1)) * lx1 + + (-sin(thetaZ1) * cos(thetaX1) + cos(thetaZ1) * sin(thetaY1) * sin(thetaX1)) * ly1 + + (sin(thetaZ1) * sin(thetaX1) + cos(thetaZ1) * sin(thetaY1) * cos(thetaX1)) * lz1; + + double cy1 = (sin(thetaZ1) * cos(thetaY1)) * lx1 + + (cos(thetaZ1) * cos(thetaX1) + sin(thetaZ1) * sin(thetaY1) * sin(thetaX1)) * ly1 + + (-cos(thetaZ1) * sin(thetaX1) + sin(thetaZ1) * sin(thetaY1) * cos(thetaX1)) * lz1; + + double cz1 = (-sin(thetaY1)) * lx1 + + (cos(thetaY1) * sin(thetaX1)) * ly1 + + (cos(thetaY1) * cos(thetaX1)) * lz1; + + double cx2 = (cos(thetaZ2) * cos(thetaY2)) * lx2 + + (-sin(thetaZ2) * cos(thetaX2) + cos(thetaZ2) * sin(thetaY2) * sin(thetaX2)) * ly2 + + (sin(thetaZ2) * sin(thetaX2) + cos(thetaZ2) * sin(thetaY2) * cos(thetaX2)) * lz2; + + double cy2 = (sin(thetaZ2) * cos(thetaY2)) * lx2 + + (cos(thetaZ2) * cos(thetaX2) + sin(thetaX2) * sin(thetaY2) * sin(thetaZ2)) * ly2 + + (-cos(thetaZ2) * sin(thetaX2) + sin(thetaZ2) * sin(thetaY2) * cos(thetaX2)) * lz2; + + double cz2 = (-sin(thetaY2)) * lx2 + + (cos(thetaY2) * sin(thetaX2)) * ly2 + + (cos(thetaY2) * cos(thetaX2)) * lz2; + //Normal displacement jump in x-direction + //First node + answer.at(1, 1) = -1.; + answer.at(1, 2) = 0.; + answer.at(1, 3) = 0.; + answer.at(1, 4) = 0.; + answer.at(1, 5) = -cz1; + answer.at(1, 6) = cy1; + //Second node + answer.at(1, 7) = 1.; + answer.at(1, 8) = 0.; + answer.at(1, 9) = 0.; + answer.at(1, 10) = 0.; + answer.at(1, 11) = -cz2; + answer.at(1, 12) = cy2; + + //Shear displacement jump in y-plane + //first node + answer.at(2, 1) = 0.; + answer.at(2, 2) = -1.; + answer.at(2, 3) = 0.; + answer.at(2, 4) = cz1; + answer.at(2, 5) = 0; + answer.at(2, 6) = -cx1; + //Second node + answer.at(2, 7) = 0.; + answer.at(2, 8) = 1.; + answer.at(2, 9) = 0.; + answer.at(2, 10) = cz2; + answer.at(2, 11) = 0; + answer.at(2, 12) = -cx2; + + //Shear displacement jump in z-plane + //first node + answer.at(3, 1) = 0.; + answer.at(3, 2) = 0.; + answer.at(3, 3) = -1.; + answer.at(3, 4) = -cy1; + answer.at(3, 5) = cx1; + answer.at(3, 6) = 0.; + //Second node + answer.at(3, 7) = 0.; + answer.at(3, 8) = 0.; + answer.at(3, 9) = 1.; + answer.at(3, 10) = -cy2; + answer.at(3, 11) = cx2; + answer.at(3, 12) = 0.; + + //Rotation around x-axis + //First node + answer.at(4, 1) = 0.; + answer.at(4, 2) = 0; + answer.at(4, 3) = 0.; + answer.at(4, 4) = -1.; + answer.at(4, 5) = 0.; + answer.at(4, 6) = 0.; + //Second node + answer.at(4, 7) = 0.; + answer.at(4, 8) = 0.; + answer.at(4, 9) = 0.; + answer.at(4, 10) = 1.; + answer.at(4, 11) = 0.; + answer.at(4, 12) = 0.; + + //Rotation around y-axis + //First node + answer.at(5, 1) = 0.; + answer.at(5, 2) = 0.; + answer.at(5, 3) = 0.; + answer.at(5, 4) = 0.; + answer.at(5, 5) = -1.; + answer.at(5, 6) = 0.; + //Second node + answer.at(5, 7) = 0.; + answer.at(5, 8) = 0.; + answer.at(5, 9) = 0.; + answer.at(5, 10) = 0.; + answer.at(5, 11) = 1.; + answer.at(5, 12) = 0.; + + //Rotation around z-axis + //First node + answer.at(6, 1) = 0.; + answer.at(6, 2) = 0.; + answer.at(6, 3) = 0.; + answer.at(6, 4) = 0.; + answer.at(6, 5) = 0.; + answer.at(6, 6) = -1.; + //Second node + answer.at(6, 7) = 0.; + answer.at(6, 8) = 0.; + answer.at(6, 9) = 0.; + answer.at(6, 10) = 0.; + answer.at(6, 11) = 0.; + answer.at(6, 12) = 1.; + + return; + } + void + LatticeFrame3dNL2::computeStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, + TimeStep *tStep) + { + FloatMatrix d, bt, db, b; + FloatArray u; + + this->computeVectorOf(VM_Total, tStep, u); + this->length = computeLength(); + + answer.resize(12, 12); + answer.zero(); + this->computeBmatrixAt(integrationRulesArray [ 0 ]->getIntegrationPoint(0), b); + + this->computeConstitutiveMatrixAt(d, rMode, integrationRulesArray [ 0 ]->getIntegrationPoint(0), tStep); + + //Rotate constitutive stiffness matrix + FloatMatrix r(6, 6), rT(6, 6), dR(6, 6), rTDR(6, 6); + computeGtoLStrainRotationMatrix(r); + rT.beTranspositionOf(r); + + dR.beProductOf(d, r); + rTDR.beProductOf(rT, dR); + + db.beProductOf(rTDR, b); + db.times(1. / length); + bt.beTranspositionOf(b); + answer.beProductOf(bt, db); + + return; + } + + void + LatticeFrame3dNL2::computeStrainVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep) + { + //Compute normalised displacement jumps in global coordinate system first and then rotate it to the local coordiante system. + + FloatArray u; + this->computeVectorOf(VM_Total, tStep, u); + + FloatArray coordA(3), coordB(3); + coordA = giveNode(1)->giveCoordinates(); + coordB = giveNode(2)->giveCoordinates(); + + this->length = computeLength(); + + double l1 = this->length * ( 1. + this->s ) / 2.; + double l2 = this->length * ( 1. - this->s ) / 2.; + + double lx1 = ( coordB.at(1) - coordA.at(1) ) * ( 1. + this->s ) / 2.; + double ly1 = ( coordB.at(2) - coordA.at(2) ) * ( 1. + this->s ) / 2.; + double lz1 = ( coordB.at(3) - coordA.at(3) ) * ( 1. + this->s ) / 2.; + + double lx2 = ( coordB.at(1) - coordA.at(1) ) * ( 1. - this->s ) / 2.; + double ly2 = ( coordB.at(2) - coordA.at(2) ) * ( 1. - this->s ) / 2.; + double lz2 = ( coordB.at(3) - coordA.at(3) ) * ( 1. - this->s ) / 2.; + + LatticeMaterialStatus *lmatStat = dynamic_cast < LatticeMaterialStatus * > ( integrationRulesArray [ 0 ]->getIntegrationPoint(0)->giveMaterialStatus() ); + auto strain = lmatStat->giveLatticeStrain(); + + //First node rotations + double uX1 = u.at(1); + double uY1 = u.at(2); + double uZ1 = u.at(3); + + double thetaX1 = u.at(4); + double thetaY1 = u.at(5); + double thetaZ1 = u.at(6); + + double thetaX2 = u.at(10); + double thetaY2 = u.at(11); + double thetaZ2 = u.at(12); + +// double cx1 = lz1 * sin(thetaY1) + lx1 * cos(thetaY1) * cos(thetaZ1) - ly1 * cos(thetaY1) * sin(thetaZ1); +// double cx2 = lz2 * sin(thetaY2) + lx2 * cos(thetaY2) * cos(thetaZ2) - ly2 * cos(thetaY2) * sin(thetaZ2); +// double cy1 = lx1 * ( cos(thetaX1) * sin(thetaZ1) + cos(thetaZ1) * sin(thetaX1) * sin(thetaY1) ) + ly1 * ( cos(thetaX1) * cos(thetaZ1) - sin(thetaX1) * sin(thetaY1) * sin(thetaZ1) ) - lz1 * cos(thetaY1) * sin(thetaX1); +// double cy2 = lx2 * ( cos(thetaX2) * sin(thetaZ2) + cos(thetaZ2) * sin(thetaX2) * sin(thetaY2) ) + ly2 * ( cos(thetaX2) * cos(thetaZ2) - sin(thetaX2) * sin(thetaY2) * sin(thetaZ2) ) - lz2 * cos(thetaY2) * sin(thetaX2); +// double cz1 = lx1 * ( sin(thetaX1) * sin(thetaZ1) - cos(thetaX1) * cos(thetaZ1) * sin(thetaY1) ) + ly1 * ( cos(thetaZ1) * sin(thetaX1) + cos(thetaX1) * sin(thetaY1) * sin(thetaZ1) ) + lz1 * cos(thetaX1) * cos(thetaY1); +// double cz2 = lx2 * ( sin(thetaX2) * sin(thetaZ2) - cos(thetaX2) * cos(thetaZ2) * sin(thetaY2) ) + ly2 * ( cos(thetaZ2) * sin(thetaX2) + cos(thetaX2) * sin(thetaY2) * sin(thetaZ2) ) + lz2 * cos(thetaX2) * cos(thetaY2); + + double cx1 = (cos(thetaZ1) * cos(thetaY1)) * lx1 + + (-sin(thetaZ1) * cos(thetaX1) + cos(thetaZ1) * sin(thetaY1) * sin(thetaX1)) * ly1 + + (sin(thetaZ1) * sin(thetaX1) + cos(thetaZ1) * sin(thetaY1) * cos(thetaX1)) * lz1; + + double cy1 = (sin(thetaZ1) * cos(thetaY1)) * lx1 + + (cos(thetaZ1) * cos(thetaX1) + sin(thetaX1) * sin(thetaY1) * sin(thetaZ1)) * ly1 + + (-cos(thetaZ1) * sin(thetaX1) + sin(thetaZ1) * sin(thetaY1) * cos(thetaX1)) * lz1; + + double cz1 = (-sin(thetaY1)) * lx1 + + (cos(thetaY1) * sin(thetaX1)) * ly1 + + (cos(thetaY1) * cos(thetaX1)) * lz1; + + double cx2 = (cos(thetaZ2) * cos(thetaY2)) * lx2 + + (-sin(thetaZ2) * cos(thetaX2) + cos(thetaZ2) * sin(thetaY2) * sin(thetaX2)) * ly2 + + (sin(thetaZ2) * sin(thetaX2) + cos(thetaZ2) * sin(thetaY2) * cos(thetaX2)) * lz2; + + double cy2 = (sin(thetaZ2) * cos(thetaY2)) * lx2 + + (cos(thetaZ2) * cos(thetaX2) + sin(thetaX2) * sin(thetaY2) * sin(thetaZ2)) * ly2 + + (-cos(thetaZ2) * sin(thetaX2) + sin(thetaZ2) * sin(thetaY2) * cos(thetaX2)) * lz2; + + double cz2 = (-sin(thetaY2)) * lx2 + + (cos(thetaY2) * sin(thetaX2)) * ly2 + + (cos(thetaY2) * cos(thetaX2)) * lz2; + + answer.resize(6); + answer.at(1) = u.at(7) - u.at(1) - cx1 - cx2 + lx1 + lx2; + answer.at(2) = u.at(8) - u.at(2) - cy1 - cy2 + ly1 + ly2; + answer.at(3) = u.at(9) - u.at(3) - cz1 - cz2 + lz1 + lz2; + answer.at(4) = u.at(10) - u.at(4); + answer.at(5) = u.at(11) - u.at(5); + answer.at(6) = u.at(12) - u.at(6); + + answer.times(1. / this->length); + + //Rotate strain vector to local coordinate system + FloatMatrix rotationMatrix(6, 6); + computeGtoLStrainRotationMatrix(rotationMatrix); + answer.rotatedWith(rotationMatrix, 'n'); + } + + + bool + LatticeFrame3dNL2::computeGtoLStrainRotationMatrix(FloatMatrix &answer) + { + FloatMatrix lcs; + answer.resize(6, 6); + answer.zero(); + + this->giveLocalCoordinateSystem(lcs); + for ( int i = 1; i <= 3; i++ ) { + for ( int j = 1; j <= 3; j++ ) { + answer.at(i, j) = lcs.at(i, j); + answer.at(i + 3, j + 3) = lcs.at(i, j); + } + } + + return 1; + } + + + bool + LatticeFrame3dNL2::computeGtoLRotationMatrix(FloatMatrix &answer) + { + return false; + } + + + + // + void + LatticeFrame3dNL2::giveInternalForcesVector(FloatArray &answer, + TimeStep *tStep, int useUpdatedGpRecord) + { + FloatMatrix b, bt, bf, d; + FloatArray u, stress, strain; + this->computeVectorOf(VM_Total, tStep, u); + this->length = computeLength(); + GaussPoint *gp = this->integrationRulesArray [ 0 ]->getIntegrationPoint(0); + + FloatArray coordA(3), coordB(3); + coordA = giveNode(1)->giveCoordinates(); + coordB = giveNode(2)->giveCoordinates(); + + this->computeStrainVector(strain, gp, tStep); + this->computeStressVector(stress, strain, integrationRulesArray [ 0 ]->getIntegrationPoint(0), tStep); + + /* // Rotate strain vector to local coordinate system */ + FloatMatrix rotationMatrix(6, 6); + computeGtoLStrainRotationMatrix(rotationMatrix); +// bt.beTranspositionOf(rotationMatrix); + stress.rotatedWith(rotationMatrix, 't'); + + //First node rotations + double thetaX1 = u.at(4); + double thetaY1 = u.at(5); + double thetaZ1 = u.at(6); + + //Second node rotations + double thetaX2 = u.at(10); + double thetaY2 = u.at(11); + double thetaZ2 = u.at(12); + + + LatticeMaterialStatus *lmatStat = dynamic_cast < LatticeMaterialStatus * > ( integrationRulesArray [ 0 ]->getIntegrationPoint(0)->giveMaterialStatus() ); + + double l1 = this->length * ( 1. + this->s ) / 2.; + double l2 = this->length * ( 1. - this->s ) / 2.; + +// + double lx1 = ( coordB.at(1) - coordA.at(1) ) * ( 1. + this->s ) / 2.; + double ly1 = ( coordB.at(2) - coordA.at(2) ) * ( 1. + this->s ) / 2.; + double lz1 = ( coordB.at(3) - coordA.at(3) ) * ( 1. + this->s ) / 2.; + + double lx2 = ( coordB.at(1) - coordA.at(1) ) * ( 1. - this->s ) / 2.; + double ly2 = ( coordB.at(2) - coordA.at(2) ) * ( 1. - this->s ) / 2.; + double lz2 = ( coordB.at(3) - coordA.at(3) ) * ( 1. - this->s ) / 2.; + + +// double cx1 = lz1 * sin(thetaY1) + lx1 * cos(thetaY1) * cos(thetaZ1) - ly1 * cos(thetaY1) * sin(thetaZ1); +// double cx2 = lz2 * sin(thetaY2) + lx2 * cos(thetaY2) * cos(thetaZ2) - ly2 * cos(thetaY2) * sin(thetaZ2); +// double cy1 = lx1 * ( cos(thetaX1) * sin(thetaZ1) + cos(thetaZ1) * sin(thetaX1) * sin(thetaY1) ) + ly1 * ( cos(thetaX1) * cos(thetaZ1) - sin(thetaX1) * sin(thetaY1) * sin(thetaZ1) ) - lz1 * cos(thetaY1) * sin(thetaX1); +// double cy2 = lx2 * ( cos(thetaX2) * sin(thetaZ2) + cos(thetaZ2) * sin(thetaX2) * sin(thetaY2) ) + ly2 * ( cos(thetaX2) * cos(thetaZ2) - sin(thetaX2) * sin(thetaY2) * sin(thetaZ2) ) - lz2 * cos(thetaY2) * sin(thetaX2); +// double cz1 = lx1 * ( sin(thetaX1) * sin(thetaZ1) - cos(thetaX1) * cos(thetaZ1) * sin(thetaY1) ) + ly1 * ( cos(thetaZ1) * sin(thetaX1) + cos(thetaX1) * sin(thetaY1) * sin(thetaZ1) ) + lz1 * cos(thetaX1) * cos(thetaY1); +// double cz2 = lx2 * ( sin(thetaX2) * sin(thetaZ2) - cos(thetaX2) * cos(thetaZ2) * sin(thetaY2) ) + ly2 * ( cos(thetaZ2) * sin(thetaX2) + cos(thetaX2) * sin(thetaY2) * sin(thetaZ2) ) + lz2 * cos(thetaX2) * cos(thetaY2); + + double cx1 = (cos(thetaZ1) * cos(thetaY1)) * lx1 + + (-sin(thetaZ1) * cos(thetaX1) + cos(thetaZ1) * sin(thetaY1) * sin(thetaX1)) * ly1 + + (sin(thetaZ1) * sin(thetaX1) + cos(thetaZ1) * sin(thetaY1) * cos(thetaX1)) * lz1; + + double cy1 = (sin(thetaZ1) * cos(thetaY1)) * lx1 + + (cos(thetaZ1) * cos(thetaX1) + sin(thetaX1) * sin(thetaY1) * sin(thetaZ1)) * ly1 + + (-cos(thetaZ1) * sin(thetaX1) + sin(thetaZ1) * sin(thetaY1) * cos(thetaX1)) * lz1; + + double cz1 = (-sin(thetaY1)) * lx1 + + (cos(thetaY1) * sin(thetaX1)) * ly1 + + (cos(thetaY1) * cos(thetaX1)) * lz1; + + double cx2 = (cos(thetaZ2) * cos(thetaY2)) * lx2 + + (-sin(thetaZ2) * cos(thetaX2) + cos(thetaZ2) * sin(thetaY2) * sin(thetaX2)) * ly2 + + (sin(thetaZ2) * sin(thetaX2) + cos(thetaZ2) * sin(thetaY2) * cos(thetaX2)) * lz2; + + double cy2 = (sin(thetaZ2) * cos(thetaY2)) * lx2 + + (cos(thetaZ2) * cos(thetaX2) + sin(thetaX2) * sin(thetaY2) * sin(thetaZ2)) * ly2 + + (-cos(thetaZ2) * sin(thetaX2) + sin(thetaZ2) * sin(thetaY2) * cos(thetaX2)) * lz2; + + double cz2 = (-sin(thetaY2)) * lx2 + + (cos(thetaY2) * sin(thetaX2)) * ly2 + + (cos(thetaY2) * cos(thetaX2)) * lz2; + + answer.resize(12); + answer.at(1) = -stress.at(1); + answer.at(2) = -stress.at(2); + answer.at(3) = -stress.at(3); + answer.at(4) = stress.at(2) * cz1 - stress.at(3) * cy1 - stress.at(4); + answer.at(5) = -stress.at(1) * cz1 + stress.at(3) * cx1 - stress.at(5); + answer.at(6) = stress.at(1) * cy1 - stress.at(2) * cx1 - stress.at(6); + answer.at(7) = stress.at(1); + answer.at(8) = stress.at(2); + answer.at(9) = stress.at(3); + answer.at(10) = stress.at(2) * cz2 - stress.at(3) * cy2 + stress.at(4); + answer.at(11) = -stress.at(1) * cz2 + stress.at(3) * cx2 + stress.at(5); + answer.at(12) = stress.at(1) * cy2 - stress.at(2) * cx2 + stress.at(6); + + + lmatStat->letTempInternalForcesBe(answer); + } +} // end namespace oofem diff --git a/src/sm/Elements/LatticeElements/latticeframe3dnl3.C b/src/sm/Elements/LatticeElements/latticeframe3dnl3.C new file mode 100644 index 000000000..d66212e9a --- /dev/null +++ b/src/sm/Elements/LatticeElements/latticeframe3dnl3.C @@ -0,0 +1,412 @@ +/* + * + * ##### ##### ###### ###### ### ### + * ## ## ## ## ## ## ## ### ## + * ## ## ## ## #### #### ## # ## + * ## ## ## ## ## ## ## ## + * ## ## ## ## ## ## ## ## + * ##### ##### ## ###### ## ## + * + * + * OOFEM : Object Oriented Finite Element Code + * + * Copyright (C) 1993 - 2023 Borek Patzak + * + * + * + * Czech Technical University, Faculty of Civil Engineering, + * Department of Structural Mechanics, 166 29 Prague, Czech Republic + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + */ + +#include "domain.h" +#include "latticeframe3dnl3.h" +#include "../sm/Materials/LatticeMaterials/latticematstatus.h" +#include "node.h" +#include "material.h" +#include "gausspoint.h" +#include "gaussintegrationrule.h" +#include "floatmatrix.h" +#include "floatmatrixf.h" +#include "intarray.h" +#include "floatarray.h" +#include "floatarrayf.h" +#include "mathfem.h" +#include "latticeframe3d.h" +#include "contextioerr.h" +#include "datastream.h" +#include "classfactory.h" +#include "sm/CrossSections/latticecrosssection.h" +#include "engngm.h" + +#ifdef __OOFEG + #include "oofeggraphiccontext.h" + #include "../sm/Materials/structuralmaterial.h" +#endif + +namespace oofem { + REGISTER_Element(LatticeFrame3dNL3); + + LatticeFrame3dNL3::LatticeFrame3dNL3(int n, Domain *aDomain) : LatticeFrame3d(n, aDomain) + { + numberOfDofMans = 2; + } + + LatticeFrame3dNL3::~LatticeFrame3dNL3() + {} + void + LatticeFrame3dNL3::computeBmatrixAt(GaussPoint *aGaussPoint, FloatMatrix &answer, int li, int ui) + //Returns the strain matrix of the receiver. + { + FloatArray u; + TimeStep *tStep = this->domain->giveEngngModel()->giveCurrentStep(); + this->computeVectorOf(VM_Total, tStep, u); + + FloatArray uGlobal(12); + FloatMatrix r(12, 12), rT(12, 12); + computeGtoLRotationMatrix(r); + rT.beTranspositionOf(r); + uGlobal.beProductOf(rT, u); + + FloatArray uROneGlobal(3); + uROneGlobal.at(1) = uGlobal.at(4); + uROneGlobal.at(2) = uGlobal.at(5); + uROneGlobal.at(3) = uGlobal.at(6); + + FloatArray uRTwoGlobal(3); + uRTwoGlobal.at(1) = uGlobal.at(10); + uRTwoGlobal.at(2) = uGlobal.at(11); + uRTwoGlobal.at(3) = uGlobal.at(12); + + //Compute the two rotation matrices in the local coordinate system + FloatMatrix rotationMatrixOne(3, 3), rotationMatrixTwo(3, 3); + ; + computeLocalRotationMatrix(rotationMatrixOne, uROneGlobal); + computeLocalRotationMatrix(rotationMatrixTwo, uRTwoGlobal); + + this->length = computeLength(); + double l1 = this->length * ( 1. + this->s ) / 2; + double l2 = this->length * ( 1. - this->s ) / 2; + LatticeMaterialStatus *lmatStat = dynamic_cast < LatticeMaterialStatus * > ( integrationRulesArray [ 0 ]->getIntegrationPoint(0)->giveMaterialStatus() ); + + double cx1 = rotationMatrixOne.at(1, 1) * l1; + double cy1 = rotationMatrixOne.at(2, 1) * l1; + double cz1 = rotationMatrixOne.at(3, 1) * l1; + + double cx2 = rotationMatrixTwo.at(1, 1) * l2; + double cy2 = rotationMatrixTwo.at(2, 1) * l2; + double cz2 = rotationMatrixTwo.at(3, 1) * l2; + + answer.resize(6, 12); + //Normal displacement jump in x-direction + //First node + answer.at(1, 1) = -1.; + answer.at(1, 2) = 0.; + answer.at(1, 3) = 0.; + answer.at(1, 4) = 0.; + answer.at(1, 5) = -cz1; + answer.at(1, 6) = cy1; + //Second node + answer.at(1, 7) = 1.; + answer.at(1, 8) = 0.; + answer.at(1, 9) = 0.; + answer.at(1, 10) = 0.; + answer.at(1, 11) = -cz2; + answer.at(1, 12) = cy2; + + //Shear displacement jump in y-plane + //first node + answer.at(2, 1) = 0.; + answer.at(2, 2) = -1.; + answer.at(2, 3) = 0.; + answer.at(2, 4) = cz1; + answer.at(2, 5) = 0; + answer.at(2, 6) = -cx1; + //Second node + answer.at(2, 7) = 0.; + answer.at(2, 8) = 1.; + answer.at(2, 9) = 0.; + answer.at(2, 10) = cz2; + answer.at(2, 11) = 0; + answer.at(2, 12) = -cx2; + + //Shear displacement jump in z-plane + //first node + answer.at(3, 1) = 0.; + answer.at(3, 2) = 0.; + answer.at(3, 3) = -1.; + answer.at(3, 4) = -cy1; + answer.at(3, 5) = cx1; + answer.at(3, 6) = 0.; + //Second node + answer.at(3, 7) = 0.; + answer.at(3, 8) = 0.; + answer.at(3, 9) = 1.; + answer.at(3, 10) = -cy2; + answer.at(3, 11) = cx2; + answer.at(3, 12) = 0.; + + //Rotation around x-axis + //First node + answer.at(4, 1) = 0.; + answer.at(4, 2) = 0; + answer.at(4, 3) = 0.; + answer.at(4, 4) = -1.; + answer.at(4, 5) = 0.; + answer.at(4, 6) = 0.; + //Second node + answer.at(4, 7) = 0.; + answer.at(4, 8) = 0.; + answer.at(4, 9) = 0.; + answer.at(4, 10) = 1.; + answer.at(4, 11) = 0.; + answer.at(4, 12) = 0.; + + //Rotation around y-axis + //First node + answer.at(5, 1) = 0.; + answer.at(5, 2) = 0.; + answer.at(5, 3) = 0.; + answer.at(5, 4) = 0.; + answer.at(5, 5) = -1.; + answer.at(5, 6) = 0.; + //Second node + answer.at(5, 7) = 0.; + answer.at(5, 8) = 0.; + answer.at(5, 9) = 0.; + answer.at(5, 10) = 0.; + answer.at(5, 11) = 1.; + answer.at(5, 12) = 0.; + + //Rotation around z-axis + //First node + answer.at(6, 1) = 0.; + answer.at(6, 2) = 0.; + answer.at(6, 3) = 0.; + answer.at(6, 4) = 0.; + answer.at(6, 5) = 0.; + answer.at(6, 6) = -1.; + //Second node + answer.at(6, 7) = 0.; + answer.at(6, 8) = 0.; + answer.at(6, 9) = 0.; + answer.at(6, 10) = 0.; + answer.at(6, 11) = 0.; + answer.at(6, 12) = 1.; + + return; + } + void + LatticeFrame3dNL3::computeStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, + TimeStep *tStep) + { + FloatMatrix d, bt, db, b; + FloatArray u; + + this->computeVectorOf(VM_Total, tStep, u); + this->length = computeLength(); + + answer.resize(12, 12); + answer.zero(); + this->computeBmatrixAt(integrationRulesArray [ 0 ]->getIntegrationPoint(0), b); + this->computeConstitutiveMatrixAt(d, rMode, integrationRulesArray [ 0 ]->getIntegrationPoint(0), tStep); + + db.beProductOf(d, b); + db.times(1. / length); + bt.beTranspositionOf(b); + answer.beProductOf(bt, db); + + return; + } + + void LatticeFrame3dNL3::computeLocalRotationMatrix(FloatMatrix &answer, FloatArray &rotation) { + //Use global rotational DOFs to compute rotation matrix. + //The order for the global rotations is first z, then y and finally x. + double thetaX = rotation.at(1); + double thetaY = rotation.at(2); + double thetaZ = rotation.at(3);; + +// +// FloatMatrix globalR(3, 3); +// globalR.at(1, 1) = cos(thetaY) * cos(thetaZ); +// globalR.at(1, 2) = -cos(thetaY) * sin(thetaZ); +// globalR.at(1, 3) = sin(thetaY); +// globalR.at(2, 1) = cos(thetaX) * sin(thetaZ) + sin(thetaX) * sin(thetaY) * cos(thetaZ); +// globalR.at(2, 2) = -sin(thetaX) * sin(thetaY) * sin(thetaZ) + cos(thetaX) * cos(thetaZ); +// globalR.at(2, 3) = -sin(thetaX) * cos(thetaY); +// globalR.at(3, 1) = sin(thetaX) * sin(thetaZ) - cos(thetaX) * sin(thetaY) * cos(thetaZ); +// globalR.at(3, 2) = cos(thetaX) * sin(thetaY) * sin(thetaZ) + sin(thetaX) * cos(thetaZ); +// globalR.at(3, 3) = cos(thetaX) * cos(thetaY); + + FloatMatrix globalR(3, 3); + globalR.at(1, 1) = cos(thetaZ) * cos(thetaY); + globalR.at(1, 2) = - sin(thetaZ)*cos(thetaX) + cos(thetaZ)*sin(thetaY)*sin(thetaX); + globalR.at(1, 3) = sin(thetaZ) * sin(thetaX) + cos(thetaZ) * sin(thetaY) * cos(thetaX); + globalR.at(2, 1) = sin(thetaZ) * cos(thetaY); + globalR.at(2, 2) = cos(thetaZ) * cos(thetaX) + sin(thetaX) * sin(thetaY) * sin(thetaZ); + globalR.at(2, 3) = - cos(thetaZ) * sin(thetaX)+ sin(thetaZ) * sin(thetaY) * cos(thetaX) ; + globalR.at(3, 1) = -sin(thetaY); + globalR.at(3, 2) = cos(thetaY) * sin(thetaX); + globalR.at(3, 3) = cos(thetaY) * cos(thetaX); + + FloatMatrix transform(3, 3), transformT(3, 3), localR(3, 3), help(3, 3); + this->giveLocalCoordinateSystem(transform); + transformT.beTranspositionOf(transform); + + help.beProductOf(globalR, transformT); + answer.beProductOf(transform, help); + + return; + } + + void + LatticeFrame3dNL3::computeStrainVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep) + { + FloatArray u; + this->computeVectorOf(VM_Total, tStep, u); + + FloatArray uGlobal(12); + FloatMatrix r(12, 12), rT(12, 12); + computeGtoLRotationMatrix(r); + rT.beTranspositionOf(r); + uGlobal.beProductOf(rT, u); + + FloatArray uROneGlobal(3); + uROneGlobal.at(1) = uGlobal.at(4); + uROneGlobal.at(2) = uGlobal.at(5); + uROneGlobal.at(3) = uGlobal.at(6); + + FloatArray uRTwoGlobal(3); + uRTwoGlobal.at(1) = uGlobal.at(10); + uRTwoGlobal.at(2) = uGlobal.at(11); + uRTwoGlobal.at(3) = uGlobal.at(12); + + //Compute the two rotation matrices in the local coordinate system + FloatMatrix rotationMatrixOne(3, 3), rotationMatrixTwo(3, 3); + + computeLocalRotationMatrix(rotationMatrixOne, uROneGlobal); + computeLocalRotationMatrix(rotationMatrixTwo, uRTwoGlobal); + + this->length = computeLength(); + double l1 = this->length * ( 1. + this->s ) / 2; + double l2 = this->length * ( 1. - this->s ) / 2; + LatticeMaterialStatus *lmatStat = dynamic_cast < LatticeMaterialStatus * > ( integrationRulesArray [ 0 ]->getIntegrationPoint(0)->giveMaterialStatus() ); + + double cx1 = rotationMatrixOne.at(1, 1) * l1; + double cy1 = rotationMatrixOne.at(2, 1) * l1; + double cz1 = rotationMatrixOne.at(3, 1) * l1; + + double cx2 = rotationMatrixTwo.at(1, 1) * l2; + double cy2 = rotationMatrixTwo.at(2, 1) * l2; + double cz2 = rotationMatrixTwo.at(3, 1) * l2; + + // + answer.resize(6); + answer.at(1) = u.at(7) - u.at(1) - cx2 - cx1 + l1 + l2; + answer.at(2) = u.at(8) - u.at(2) - cy2 - cy1; + answer.at(3) = u.at(9) - u.at(3) - cz2 - cz1; + answer.at(4) = u.at(10) - u.at(4); + answer.at(5) = u.at(11) - u.at(5); + answer.at(6) = u.at(12) - u.at(6); + answer.times(1. / this->length); + } + + + bool + LatticeFrame3dNL3::computeGtoLComponentTransformationMatrix(FloatMatrix &answer) + { + FloatMatrix lcs; + answer.resize(3, 3); + answer.zero(); + + this->giveLocalCoordinateSystem(lcs); + for ( int i = 1; i <= 3; i++ ) { + for ( int j = 1; j <= 3; j++ ) { + answer.at(i, j) = lcs.at(i, j); + } + } + + return 1; + } + + + void + LatticeFrame3dNL3::giveInternalForcesVector(FloatArray &answer, + TimeStep *tStep, int useUpdatedGpRecord) + { + FloatMatrix b, bt, bf, d; + FloatArray u, stress, strain; + this->computeVectorOf(VM_Total, tStep, u); + + GaussPoint *gp = this->integrationRulesArray [ 0 ]->getIntegrationPoint(0); + + this->computeStrainVector(strain, gp, tStep); + this->computeStressVector(stress, strain, integrationRulesArray [ 0 ]->getIntegrationPoint(0), tStep); + + LatticeMaterialStatus *lmatStat = dynamic_cast < LatticeMaterialStatus * > ( integrationRulesArray [ 0 ]->getIntegrationPoint(0)->giveMaterialStatus() ); + + + FloatArray uGlobal(12); + FloatMatrix r(12, 12), rT(12, 12); + computeGtoLRotationMatrix(r); + rT.beTranspositionOf(r); + uGlobal.beProductOf(rT, u); + + FloatArray uROneGlobal(3); + uROneGlobal.at(1) = uGlobal.at(4); + uROneGlobal.at(2) = uGlobal.at(5); + uROneGlobal.at(3) = uGlobal.at(6); + + FloatArray uRTwoGlobal(3); + uRTwoGlobal.at(1) = uGlobal.at(10); + uRTwoGlobal.at(2) = uGlobal.at(11); + uRTwoGlobal.at(3) = uGlobal.at(12); + + //Compute the two rotation matrices in the local coordinate system + FloatMatrix rotationMatrixOne(3, 3), rotationMatrixTwo(3, 3); + ; + computeLocalRotationMatrix(rotationMatrixOne, uROneGlobal); + computeLocalRotationMatrix(rotationMatrixTwo, uRTwoGlobal); + + this->length = computeLength(); + double l1 = this->length * ( 1. + this->s ) / 2; + double l2 = this->length * ( 1. - this->s ) / 2; + + double cx1 = rotationMatrixOne.at(1, 1) * l1; + double cy1 = rotationMatrixOne.at(2, 1) * l1; + double cz1 = rotationMatrixOne.at(3, 1) * l1; + + double cx2 = rotationMatrixTwo.at(1, 1) * l2; + double cy2 = rotationMatrixTwo.at(2, 1) * l2; + double cz2 = rotationMatrixTwo.at(3, 1) * l2; + + answer.resize(12); + answer.at(1) = -stress.at(1); + answer.at(2) = -stress.at(2); + answer.at(3) = -stress.at(3); + answer.at(4) = stress.at(2) * cz1 - stress.at(3) * cy1 - stress.at(4); + answer.at(5) = -stress.at(1) * cz1 + stress.at(3) * cx1 - stress.at(5); + answer.at(6) = stress.at(1) * cy1 - stress.at(2) * cx1 - stress.at(6); + answer.at(7) = stress.at(1); + answer.at(8) = stress.at(2); + answer.at(9) = stress.at(3); + answer.at(10) = stress.at(2) * cz2 - stress.at(3) * cy2 + stress.at(4); + answer.at(11) = -stress.at(1) * cz2 + stress.at(3) * cx2 + stress.at(5); + answer.at(12) = stress.at(1) * cy2 - stress.at(2) * cx2 + stress.at(6); + + lmatStat->letTempInternalForcesBe(answer); + } +} // end namespace oofem diff --git a/src/sm/Materials/LatticeMaterials/CMakeLists.txt b/src/sm/Materials/LatticeMaterials/CMakeLists.txt new file mode 100644 index 000000000..e69de29bb diff --git a/src/sm/Materials/LatticeMaterials/latticeframeconcreteplastic.C b/src/sm/Materials/LatticeMaterials/latticeframeconcreteplastic.C new file mode 100644 index 000000000..6d394b898 --- /dev/null +++ b/src/sm/Materials/LatticeMaterials/latticeframeconcreteplastic.C @@ -0,0 +1,934 @@ +/* + * + * ##### ##### ###### ###### ### ### + * ## ## ## ## ## ## ## ### ## + * ## ## ## ## #### #### ## # ## + * ## ## ## ## ## ## ## ## + * ## ## ## ## ## ## ## ## + * ##### ##### ## ###### ## ## + * + * + * OOFEM : Object Oriented Finite Element Code + * + * Copyright (C) 1993 - 2022 Borek Patzak + * + * + * + * + * Czech Technical University, Faculty of Civil Engineering, + * Department of Structural Mechanics, 166 29 Prague, Czech Republic + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + */ + +#include "latticeframeconcreteplastic.h" +#include "latticeframeelastic.h" +#include "latticestructuralmaterial.h" +#include "gausspoint.h" +#include "floatmatrix.h" +#include "floatmatrixf.h" +#include "floatarray.h" +#include "floatarrayf.h" +#include "CrossSections/structuralcrosssection.h" +#include "engngm.h" +#include "mathfem.h" +#include "Elements/LatticeElements/latticestructuralelement.h" +#include "datastream.h" +#include "staggeredproblem.h" +#include "contextioerr.h" +#include "classfactory.h" + +namespace oofem { +REGISTER_Material(LatticeFrameConcretePlastic); + +bool LatticeFrameConcretePlastic::hasMaterialModeCapability(MaterialMode mode) const +{ + return ( mode == _3dLattice ); +} + + +void LatticeFrameConcretePlastic::initializeFrom(InputRecord &ir) +{ + LatticeStructuralMaterial::initializeFrom(ir); + + // Young's modulus of the material that the beam element is made of + IR_GIVE_FIELD(ir, this->e, _IFT_LatticeFrameConcretePlastic_e); // Macro + + // Poisson's ratio of the material that the beam element is made of + IR_GIVE_FIELD(ir, this->nu, _IFT_LatticeFrameConcretePlastic_n); // Macro + + // Nx0 + IR_GIVE_FIELD(ir, this->nx0, _IFT_LatticeFrameConcretePlastic_nx0); // Macro + + // Nx01 + IR_GIVE_FIELD(ir, this->nx01, _IFT_LatticeFrameConcretePlastic_nx01); // Macro + + // vy0 + IR_GIVE_FIELD(ir, this->vy0, _IFT_LatticeFrameConcretePlastic_vy0); // Macro + + // vy01 + IR_GIVE_FIELD(ir, this->vy01, _IFT_LatticeFrameConcretePlastic_vy01); // Macro + + // vz0 + IR_GIVE_FIELD(ir, this->vz0, _IFT_LatticeFrameConcretePlastic_vz0); // Macro + + // vz01 + IR_GIVE_FIELD(ir, this->vz01, _IFT_LatticeFrameConcretePlastic_vz01); // Macro + + // Mx0 + IR_GIVE_FIELD(ir, this->mx0, _IFT_LatticeFrameConcretePlastic_mx0); // Macro + + // Mx01 + IR_GIVE_FIELD(ir, this->mx01, _IFT_LatticeFrameConcretePlastic_mx01); // Macro + + // My0 + IR_GIVE_FIELD(ir, this->my0, _IFT_LatticeFrameConcretePlastic_my0); // Macro + + // My01 + IR_GIVE_FIELD(ir, this->my01, _IFT_LatticeFrameConcretePlastic_my01); // Macro + + // Mz0 + IR_GIVE_FIELD(ir, this->mz0, _IFT_LatticeFrameConcretePlastic_mz0); // Macro + + // Mz01 + IR_GIVE_FIELD(ir, this->mz01, _IFT_LatticeFrameConcretePlastic_mz01); // Macro + + yieldTol = 1.e-6; + IR_GIVE_FIELD(ir, this->yieldTol, _IFT_LatticeFrameConcretePlastic_tol); // Macro + + this->newtonIter = 100; + IR_GIVE_FIELD(ir, this->newtonIter, _IFT_LatticeFrameConcretePlastic_iter); // Macro + + numberOfSubIncrements = 10; + IR_GIVE_FIELD(ir, this->numberOfSubIncrements, _IFT_LatticeFrameConcretePlastic_sub); // Macro + + this->plasticFlag = 1; + IR_GIVE_OPTIONAL_FIELD(ir, plasticFlag, _IFT_LatticeFrameConcretePlastic_plastic); // Macro + + // wu + IR_GIVE_OPTIONAL_FIELD(ir, this->wu, _IFT_LatticeFrameConcretePlastic_wu); + + // wf + IR_GIVE_FIELD(ir, wf, _IFT_LatticeFrameConcretePlastic_wf); +} +MaterialStatus * +LatticeFrameConcretePlastic::CreateStatus(GaussPoint *gp) const +{ + return new LatticeFrameConcretePlasticStatus(1, LatticeFrameConcretePlastic::domain, gp); +} + +MaterialStatus * +LatticeFrameConcretePlastic::giveStatus(GaussPoint *gp) const +{ + MaterialStatus *status = static_cast< MaterialStatus * >( gp->giveMaterialStatus() ); + if ( !status ) { + // create a new one + status = this->CreateStatus(gp); + + if ( status ) { + gp->setMaterialStatus(status); + } + } + // + return status; +} +int +LatticeFrameConcretePlastic::giveIPValue(FloatArray &answer, + GaussPoint *gp, + InternalStateType type, + TimeStep *atTime) +{ + auto status = static_cast< LatticeFrameConcretePlasticStatus * >( this->giveStatus(gp) ); + + if ( type == IST_PlasticLatticeStrain ) { + answer.resize(6); + answer.zero(); + answer = status->givePlasticLatticeStrain(); + return 1; + } else if ( type == IST_DamageScalar ) { + answer.resize(1); + answer.zero(); + answer.at(1) = status->giveDamage(); + return 1; + } else { + return LatticeStructuralMaterial::giveIPValue(answer, gp, type, atTime); + } +} + +double +LatticeFrameConcretePlastic::computeYieldValue(const FloatArrayF< 6 > &stress, const FloatArrayF< 6 > &k, + GaussPoint *gp, + TimeStep *tStep) const +{ + double yieldValue = 0.; + double nx = stress.at(1); + double vy = stress.at(2); + double vz = stress.at(3); + double mx = stress.at(4); + double my = stress.at(5); + double mz = stress.at(6); + + double ax; + if ( k.at(1) == 1 ) { + ax = nx0; + } else { + ax = nx01; + } + double ay; + if ( k.at(2) == 1 ) { + ay = vy0; + } else { + ay = vy01; + } + double az; + if ( k.at(3) == 1 ) { + az = vz0; + } else { + az = vz01; + } + double bx; + + if ( k.at(4) == 1 ) { + bx = mx0; + } else { + bx = mx01; + } + double by; + + if ( k.at(5) == 1 ) { + by = my0; + } else { + by = my01; + } + + double bz; + + if ( k.at(6) == 1 ) { + bz = mz0; + } else { + bz = mz01; + } + + + + yieldValue = pow(nx / ax, 2.) + pow(vy / ay, 2.) + pow(vz / az, 2.) + pow(mx / bx, 2.) + pow(my / by, 2.) + pow(mz / bz, 2.) - 1.; + + + return yieldValue; +} + +FloatArrayF< 6 > +LatticeFrameConcretePlastic::computeFVector(const FloatArrayF< 6 > &stress, const FloatArrayF< 6 > &k, + GaussPoint *gp, + TimeStep *tStep) const +{ + double nx = stress.at(1); + double vy = stress.at(2); + double vz = stress.at(3); + double mx = stress.at(4); + double my = stress.at(5); + double mz = stress.at(6); + + double ax = 0.; + if ( k.at(1) == 1 ) { + ax = nx0; + } else { + ax = nx01; + } + + double ay = 0.; + if ( k.at(2) == 1 ) { + ay = vy0; + } else { + ay = vy01; + } + + double az = 0.; + if ( k.at(3) == 1 ) { + az = vz0; + } else { + az = vz01; + } + + double bx = 0.; + if ( k.at(4) == 1 ) { + bx = mx0; + } else { + bx = mx01; + } + + double by = 0.; + if ( k.at(5) == 1 ) { + by = my0; + } else { + by = my01; + } + + double bz = 0.; + if ( k.at(6) == 1 ) { + bz = mz0; + } else { + bz = mz01; + } + + FloatArrayF< 6 >f; + f.at(1) = 2. * nx / pow(ax, 2.); + f.at(2) = 2. * vy / pow(ay, 2.); + f.at(3) = 2. * vz / pow(az, 2.); + f.at(4) = 2. * mx / pow(bx, 2.); + f.at(5) = 2. * my / pow(by, 2.); + f.at(6) = 2. * mz / pow(bz, 2.); + + return f; +} + +FloatMatrixF< 6, 6 > +LatticeFrameConcretePlastic::computeDMMatrix(const FloatArrayF< 6 > &stress, const FloatArrayF< 6 > &k, GaussPoint *gp, TimeStep *tStep) const +{ + double ax = nx01; + if ( k.at(1) == 1 ) { + ax = nx0; + } + + double ay = vy01; + ; + if ( k.at(2) == 1 ) { + ay = vy0; + } + + double az = vz01; + if ( k.at(3) == 1 ) { + az = vz0; + } + + double bx = mx01; + if ( k.at(4) == 1 ) { + bx = mx0; + } + + double by = my01; + if ( k.at(5) == 1 ) { + by = my0; + } + + double bz = mz01; + ; + if ( k.at(6) == 1 ) { + bz = mz0; + } + + FloatArrayF< 6 >dm = { + 2. / pow(ax, 2.), + 2. / pow(ay, 2.), + 2. / pow(az, 2.), + 2. / pow(bx, 2.), + 2. / pow(by, 2.), + 2. / pow(bz, 2.) + }; + + return diag(dm); +} + +FloatArrayF< 6 > +LatticeFrameConcretePlastic::giveThermalDilatationVector(GaussPoint *gp, TimeStep *tStep) const +// returns a FloatArray(6) of initial strain vector caused by unit temperature in direction of gp (element) local axes +{ + double alpha = this->give(tAlpha, gp); + + return { + alpha, 0., 0., 0., 0., 0. + }; +} + +FloatArrayF< 6 > +LatticeFrameConcretePlastic::giveStrain(GaussPoint *gp, TimeStep *tStep) const +{ + auto status = static_cast< LatticeMaterialStatus * >( this->giveStatus(gp) ); + return status->giveLatticeStrain(); +} + +const FloatArrayF< 6 > +LatticeFrameConcretePlastic::checkTransition(const FloatArrayF< 6 > &stress, GaussPoint *gp, TimeStep *tStep) const +{ + FloatArrayF< 6 >k; + + IntArray m(6); + + if ( stress.at(1) >= 0 ) { + m.at(1) = 1; + } else { + m.at(1) = -1; + } + + if ( stress.at(2) >= 0 ) { + m.at(2) = 1; + } else { + m.at(2) = -1; + } + + if ( stress.at(3) >= 0 ) { + m.at(3) = 1; + } else { + m.at(3) = -1; + } + + if ( stress.at(4) >= 0 ) { + m.at(4) = 1; + } else { + m.at(4) = -1; + } + + if ( stress.at(5) >= 0 ) { + m.at(5) = 1; + } else { + m.at(5) = -1; + } + + if ( stress.at(6) >= 0 ) { + m.at(6) = 1; + } else { + m.at(6) = -1; + } + + int counter = 1; + for ( int xi = -1; xi < 2; xi += 2 ) { + for ( int yi = -1; yi < 2; yi += 2 ) { + for ( int zi = -1; zi < 2; zi += 2 ) { + for ( int xj = -1; xj < 2; xj += 2 ) { + for ( int yj = -1; yj < 2; yj += 2 ) { + for ( int zj = -1; zj < 2; zj += 2 ) { + if ( !( zj == 0 && yj == 0 && xj == 0 && zi == 0 && yi == 0 && xi == 0 ) ) { + if ( xi == m.at(1) && yi == m.at(2) && zi == m.at(3) && xj == m.at(4) && yj == m.at(5) && zj == m.at(6) ) { + k.at(1) = xi; + k.at(2) = yi; + k.at(3) = zi; + k.at(4) = xj; + k.at(5) = yj; + k.at(6) = zj; + } + counter++; + } + } + } + } + } + } + } + if + ( k.at(1) != m.at(1) && k.at(2) != m.at(2) && k.at(2) != m.at(2) && k.at(2) != m.at(2) && k.at(2) != m.at(2) && k.at(2) != m.at(2) ) { + OOFEM_ERROR("This case should not exist"); + } + + return k; +} + +// + +FloatArrayF< 6 > +LatticeFrameConcretePlastic::performPlasticityReturn(GaussPoint *gp, const FloatArrayF< 6 > &strain, TimeStep *tStep) const +{ + LatticeFrameConcretePlastic_ReturnResult returnResult = RR_Unknown; + int kIter = 0; + double g = this->e / ( 2. * ( 1. + this->nu ) ); + const double area = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveArea(); + const double shearareay = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveShearAreaY(); + const double shearareaz = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveShearAreaZ(); + const double ik = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveIk(); + const double iy = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveIy(); + const double iz = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveIz(); + + auto status = static_cast< LatticeFrameConcretePlasticStatus * >( this->giveStatus(gp) ); + + + /* Get plastic strain vector from status*/ + auto tempPlasticStrain = status->givePlasticLatticeStrain(); + + FloatArrayF< 6 >tangent = { area *this->e, g *shearareay, g *shearareaz, ik *g, iy *this->e, iz *this->e }; + + /* Compute trial stress*/ + auto stress = mult(tangent, strain - tempPlasticStrain); + + + + auto oldStrain = this->giveStrain(gp, tStep); + //test + auto k = checkTransition(stress, gp, tStep); + + /* Compute yield value*/ + double yieldValue = computeYieldValue(stress, k, gp, tStep); + int subIncrementCounter = 0; + + /* Check yield condition, i.e. if the yield value is less than the yield tolerance. If yield condition is valid. Do perform regular return (closest point return)*/ + + if ( yieldValue > yieldTol ) { + int subIncrementFlag = 0; + auto convergedStrain = oldStrain; + auto tempStrain = strain; + auto deltaStrain = strain - oldStrain; + // To get into the loop + status->letTempReturnResultBe(LatticeFrameConcretePlastic::RR_NotConverged); + while ( status->giveTempReturnResult() == RR_NotConverged || subIncrementFlag == 1 ) { + stress = mult(tangent, tempStrain - tempPlasticStrain); + performRegularReturn(stress, returnResult, k, yieldValue, gp, tStep); + //if ( status->giveTempReturnResult() == RR_WrongSurface ) { + // int RestartWithNewKValue; + //} + if ( status->giveTempReturnResult() == RR_NotConverged ) { + subIncrementCounter++; + if ( subIncrementCounter > numberOfSubIncrements ) { + OOFEM_LOG_INFO( "Unstable element %d \n", gp->giveElement()->giveGlobalNumber() ); + OOFEM_LOG_INFO("Yield value %e \n", yieldValue); + OOFEM_LOG_INFO( "ConvergedStrain value %e %e %e %e\n", convergedStrain.at(1), convergedStrain.at(2), convergedStrain.at(3), convergedStrain.at(4), convergedStrain.at(5), convergedStrain.at(6) ); + OOFEM_LOG_INFO( "tempStrain value %e %e %e %e\n", tempStrain.at(1), tempStrain.at(2), tempStrain.at(3), tempStrain.at(4), tempStrain.at(5), tempStrain.at(6) ); + OOFEM_LOG_INFO( "deltaStrain value %e %e %e %e\n", deltaStrain.at(1), deltaStrain.at(2), deltaStrain.at(3), deltaStrain.at(4), deltaStrain.at(5), deltaStrain.at(6) ); + OOFEM_LOG_INFO( "targetstrain value %e %e %e %e\n", strain.at(1), strain.at(2), strain.at(3), strain.at(4), strain.at(5), strain.at(6) ); + + OOFEM_ERROR("LatticeFrameConcretePlastic :: performPlasticityReturn - Could not reach convergence with small deltaStrain, giving up."); + } + subIncrementFlag = 1; + deltaStrain *= 0.5; + tempStrain = convergedStrain + deltaStrain; + } else if ( status->giveTempReturnResult() == RR_Converged && subIncrementFlag == 1 ) { + tempPlasticStrain.at(1) = tempStrain.at(1) - stress.at(1) / ( area * this->e ); + tempPlasticStrain.at(2) = tempStrain.at(2) - stress.at(2) / ( g * shearareay ); + tempPlasticStrain.at(3) = tempStrain.at(3) - stress.at(3) / ( g * shearareaz ); + tempPlasticStrain.at(4) = tempStrain.at(4) - stress.at(4) / ( ik * g ); + tempPlasticStrain.at(5) = tempStrain.at(5) - stress.at(5) / ( iy * this->e ); + tempPlasticStrain.at(6) = tempStrain.at(6) - stress.at(6) / ( iz * this->e ); + + status->letTempPlasticLatticeStrainBe( assemble< 6 >(tempPlasticStrain, { 0, 1, 2, 3, 4, 5 }) ); + + subIncrementFlag = 0; + + status->letTempReturnResultBe(LatticeFrameConcretePlastic::RR_NotConverged); + convergedStrain = tempStrain; + deltaStrain = strain - convergedStrain; + tempStrain = strain; + subIncrementCounter = 0; + } else { //Converged + //Check the surface + FloatArrayF< 6 >kCheck = checkTransition(stress, gp, tStep); + + if ( kCheck.at(1) == k.at(1) && kCheck.at(2) == k.at(2) && kCheck.at(3) == k.at(3) && kCheck.at(4) == k.at(4) && kCheck.at(5) == k.at(5) && kCheck.at(6) == k.at(6) ) { + returnResult = RR_Converged; + } else { + //ended up in a region for which other surface should have been used. + //Try with other surface + if ( kIter == 1 ) { + OOFEM_ERROR("LatticePlasticityDamage :: performPlasticityReturn - Tried both surfaces"); + } + + k = kCheck; + returnResult = RR_NotConverged; + kIter++; + } + } + } + } + + tempPlasticStrain.at(1) = strain.at(1) - stress.at(1) / ( area * this->e ); + tempPlasticStrain.at(2) = strain.at(2) - stress.at(2) / ( g * shearareay ); + tempPlasticStrain.at(3) = strain.at(3) - stress.at(3) / ( g * shearareaz ); + tempPlasticStrain.at(4) = strain.at(4) - stress.at(4) / ( ik * g ); + tempPlasticStrain.at(5) = strain.at(5) - stress.at(5) / ( iy * this->e ); + tempPlasticStrain.at(6) = strain.at(6) - stress.at(6) / ( iz * this->e ); + + status->letTempPlasticLatticeStrainBe(tempPlasticStrain); + + return stress; +} + +Interface * +LatticeFrameConcretePlastic::giveInterface(InterfaceType type) +{ + return nullptr; +} +// +void +LatticeFrameConcretePlastic::performRegularReturn(FloatArrayF< 6 > &stress, LatticeFrameConcretePlastic_ReturnResult, const FloatArrayF< 6 > &k, + double yieldValue, + GaussPoint *gp, + TimeStep *tStep) const +{ + // Use material specific status + auto status = static_cast< LatticeFrameConcretePlasticStatus * >( this->giveStatus(gp) ); + + double deltaLambda = 0.; + + auto trialStress = stress; + auto tempStress = trialStress; + + // initialise unknowns + FloatArrayF< 7 >unknowns; + unknowns.at(1) = trialStress.at(1); + unknowns.at(2) = trialStress.at(2); + unknowns.at(3) = trialStress.at(3); + unknowns.at(4) = trialStress.at(4); + unknowns.at(5) = trialStress.at(5); + unknowns.at(6) = trialStress.at(6); + unknowns.at(7) = 0.; + + yieldValue = computeYieldValue(tempStress, k, gp, tStep); + + // initiate residuals + FloatArrayF< 7 >residuals; + residuals.at(7) = yieldValue; + double normOfResiduals = 1.; // just to get into the loop + int iterationCount = 0; + while ( normOfResiduals > yieldTol ) { + iterationCount++; + if ( iterationCount == newtonIter ) { + status->letTempReturnResultBe(LatticeFrameConcretePlasticStatus::RR_NotConverged); + return; + } + + + double ax; + if ( k.at(1) == 1 ) { + ax = nx0; + } else { + ax = nx01; + } + double ay; + if ( k.at(2) == 1 ) { + ay = vy0; + } else { + ay = vy01; + } + double az; + if ( k.at(3) == 1 ) { + az = vz0; + } else { + az = vz01; + } + double bx; + + if ( k.at(4) == 1 ) { + bx = mx0; + } else { + bx = mx01; + } + double by; + + if ( k.at(5) == 1 ) { + by = my0; + } else { + by = my01; + } + + double bz; + + if ( k.at(6) == 1 ) { + bz = mz0; + } else { + bz = mz01; + } + + FloatArrayF< 7 >residualsNorm; + residualsNorm.at(1) = residuals.at(1) / ax; + residualsNorm.at(2) = residuals.at(2) / ay; + residualsNorm.at(3) = residuals.at(3) / az; + residualsNorm.at(4) = residuals.at(4) / bx; + residualsNorm.at(5) = residuals.at(5) / by; + residualsNorm.at(6) = residuals.at(6) / bz; + residualsNorm.at(7) = residuals.at(7); + normOfResiduals = norm(residualsNorm); + // printf( "normOfResiduals=%e\n", normOfResiduals ); + if ( std::isnan(normOfResiduals) ) { + status->letTempReturnResultBe(LatticeFrameConcretePlasticStatus::RR_NotConverged); + return; + } + + if ( normOfResiduals > yieldTol ) { + auto jacobian = computeJacobian(tempStress, k, deltaLambda, gp, tStep); + + auto solution = solve_check(jacobian, residuals); + if ( solution.first ) { + unknowns -= solution.second; + } else { + status->letTempReturnResultBe(LatticeFrameConcretePlastic::RR_NotConverged); + } + unknowns.at(7) = max(unknowns.at(7), 0.); // Keep deltaLambda greater than zero! + + /* Update increments final values and DeltaLambda*/ + tempStress.at(1) = unknowns.at(1); + tempStress.at(2) = unknowns.at(2); + tempStress.at(3) = unknowns.at(3); + tempStress.at(4) = unknowns.at(4); + tempStress.at(5) = unknowns.at(5); + tempStress.at(6) = unknowns.at(6); + deltaLambda = unknowns.at(7); + + /* Compute the fVector*/ + auto FVector = computeFVector(tempStress, k, gp, tStep); + double g = this->e / ( 2. * ( 1. + this->nu ) ); + const double area = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveArea(); + const double shearareay = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveShearAreaY(); + const double shearareaz = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveShearAreaZ(); + const double ik = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveIk(); + const double iy = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveIy(); + const double iz = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveIz(); + + residuals.at(1) = tempStress.at(1) - trialStress.at(1) + area * this->e * deltaLambda * FVector.at(1); + residuals.at(2) = tempStress.at(2) - trialStress.at(2) + shearareay * g * deltaLambda * FVector.at(2); + residuals.at(3) = tempStress.at(3) - trialStress.at(3) + shearareaz * g * deltaLambda * FVector.at(3); + residuals.at(4) = tempStress.at(4) - trialStress.at(4) + ik * g * deltaLambda * FVector.at(4); + residuals.at(5) = tempStress.at(5) - trialStress.at(5) + iy * this->e * deltaLambda * FVector.at(5); + residuals.at(6) = tempStress.at(6) - trialStress.at(6) + iz * this->e * deltaLambda * FVector.at(6); + residuals.at(7) = computeYieldValue(tempStress, k, gp, tStep); + } + } + status->letTempReturnResultBe(LatticeFrameConcretePlastic::RR_Converged); + stress = tempStress; +} + +FloatMatrixF< 7, 7 > +LatticeFrameConcretePlastic::computeJacobian(const FloatArrayF< 6 > &stress, const FloatArrayF< 6 > &k, + const double deltaLambda, + GaussPoint *gp, + TimeStep *tStep) const +{ + auto dMMatrix = computeDMMatrix(stress, k, gp, tStep); + auto fVector = computeFVector(stress, k, gp, tStep); + double g = this->e / ( 2. * ( 1. + this->nu ) ); + const double area = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveArea(); + const double shearareay = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveShearAreaY(); + const double shearareaz = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveShearAreaZ(); + const double ik = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveIk(); + const double iy = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveIy(); + const double iz = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveIz(); + + /* Compute matrix*/ + FloatMatrixF< 7, 7 >jacobian; + jacobian.at(1, 1) = 1. + this->e * area * deltaLambda * dMMatrix.at(1, 1); + jacobian.at(1, 2) = 0.; + jacobian.at(1, 3) = 0.; + jacobian.at(1, 4) = 0.; + jacobian.at(1, 5) = 0.; + jacobian.at(1, 6) = 0.; + jacobian.at(1, 7) = this->e * area * fVector.at(1); + + jacobian.at(2, 1) = 0.; + jacobian.at(2, 2) = 1. + shearareay * g * deltaLambda * dMMatrix.at(2, 2); + jacobian.at(2, 3) = 0.; + jacobian.at(2, 4) = 0.; + jacobian.at(2, 5) = 0.; + jacobian.at(2, 6) = 0.; + jacobian.at(2, 7) = shearareay * g * fVector.at(2); + + jacobian.at(3, 1) = 0.; + jacobian.at(3, 2) = 0.; + jacobian.at(3, 3) = 1. + shearareaz * g * deltaLambda * dMMatrix.at(3, 3); + jacobian.at(3, 4) = 0.; + jacobian.at(3, 5) = 0.; + jacobian.at(3, 6) = 0.; + jacobian.at(3, 7) = shearareaz * g * fVector.at(3); + + jacobian.at(4, 1) = 0.; + jacobian.at(4, 2) = 0.; + jacobian.at(4, 3) = 0.; + jacobian.at(4, 4) = 1. + ik * g * deltaLambda * dMMatrix.at(4, 4); + jacobian.at(4, 5) = 0.; + jacobian.at(4, 6) = 0.; + jacobian.at(4, 7) = ik * g * fVector.at(4); + + jacobian.at(5, 1) = 0.; + jacobian.at(5, 2) = 0.; + jacobian.at(5, 3) = 0.; + jacobian.at(5, 4) = 0.; + jacobian.at(5, 5) = 1. + iy * this->e * deltaLambda * dMMatrix.at(5, 5); + jacobian.at(5, 6) = 0.; + jacobian.at(5, 7) = iy * this->e * fVector.at(5); + + + jacobian.at(6, 1) = 0.; + jacobian.at(6, 2) = 0.; + jacobian.at(6, 3) = 0.; + jacobian.at(6, 4) = 0.; + jacobian.at(6, 5) = 0.; + jacobian.at(6, 6) = 1. + this->e * iz * deltaLambda * dMMatrix.at(6, 6); + jacobian.at(6, 7) = this->e * iz * fVector.at(6); + + jacobian.at(7, 1) = fVector.at(1); + jacobian.at(7, 2) = fVector.at(2); + jacobian.at(7, 3) = fVector.at(3); + jacobian.at(7, 4) = fVector.at(4); + jacobian.at(7, 5) = fVector.at(5); + jacobian.at(7, 6) = fVector.at(6); + jacobian.at(7, 7) = 0.; + + return jacobian; +} + + +FloatArrayF< 6 > +LatticeFrameConcretePlastic::giveFrameForces3d(const FloatArrayF< 6 > &originalStrain, GaussPoint *gp, TimeStep *tStep) +{ + auto status = static_cast< LatticeFrameConcretePlasticStatus * >( this->giveStatus(gp) ); + auto strain = originalStrain; + auto thermalStrain = this->computeStressIndependentStrainVector(gp, tStep, VM_Total); + if ( thermalStrain.giveSize() ) { + strain -= FloatArrayF< 6 >(thermalStrain); + } + auto stress = this->performPlasticityReturn(gp, strain, tStep); + + auto tempPlasticStrain = status->giveTempPlasticLatticeStrain(); + + double le = static_cast< LatticeStructuralElement * >( gp->giveElement() )->giveLength(); + + //double equivalentStrain = sqrt(pow(tempPlasticStrain.at(1), 2) + pow(tempPlasticStrain.at(2), 2) + pow(tempPlasticStrain.at(3), 2) + + // pow(tempPlasticStrain.at(4), 2) + pow(tempPlasticStrain.at(5), 2.) + pow(tempPlasticStrain.at(6), 2) ) - wu / le; + double equivalentStrain = sqrt(pow(tempPlasticStrain.at(4), 2) + pow(tempPlasticStrain.at(5), 2.) + pow(tempPlasticStrain.at(6), 2) ) - wu / le; + // printf("equivalentStrain = %e, wu/le = %e, le = %e\n", equivalentStrain, wu/le, le); + + double tempKappaD = 0.0; + + if ( equivalentStrain > status->giveKappaD() ) { + tempKappaD = equivalentStrain; + } else { + tempKappaD = status->giveKappaD(); + } + + double omega = 0.0; + if ( tempKappaD <= 0. ) { + omega = 0.; + } else if ( tempKappaD > 0. ) { + omega = 1. - exp(-tempKappaD * le / ( wf ) ); + } else { + printf("Should not be here\n"); + } + + if ( omega > 1.0 ) { + omega = 1.; + } else if ( omega < 0.0 ) { + omega = 0.; + } + + stress *= ( 1. - omega ); + + status->letTempLatticeStrainBe(originalStrain); + status->letTempReducedLatticeStrainBe(strain); + status->letTempLatticeStressBe(stress); + status->setTempKappaD(tempKappaD); + status->setTempDamage(omega); + + return stress; +} + + +FloatMatrixF< 6, 6 > +LatticeFrameConcretePlastic::give3dFrameStiffnessMatrix(MatResponseMode rmode, GaussPoint *gp, TimeStep *atTime) const +{ + static_cast< LatticeFrameConcretePlasticStatus * >( this->giveStatus(gp) ); + + double g = this->e / ( 2. * ( 1. + this->nu ) ); + const double area = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveArea(); + const double shearareay = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveShearAreaY(); + const double shearareaz = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveShearAreaZ(); + const double ik = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveIk(); + const double iy = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveIy(); + const double iz = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveIz(); + + + FloatArrayF< 6 >d = { + this->e * area, + g *shearareay, + g *shearareaz, + g *ik, + this->e * iy, + this->e * iz + }; + + return diag(d); +} + +LatticeFrameConcretePlasticStatus::LatticeFrameConcretePlasticStatus(int n, Domain *d, GaussPoint *g) : + LatticeMaterialStatus(g) +{} +void +LatticeFrameConcretePlasticStatus::initTempStatus() +{ + LatticeMaterialStatus::initTempStatus(); + this->tempKappaD = this->kappaD; + this->tempDamage = this->damage; +} + +void +LatticeFrameConcretePlasticStatus::updateYourself(TimeStep *atTime) +// +// updates variables (nonTemp variables describing situation at previous equilibrium state) +// after a new equilibrium state has been reached +// temporary variables are having values corresponding to newly reached equilibrium. +// +{ + LatticeMaterialStatus::updateYourself(atTime); + this->kappaD = this->tempKappaD; + this->damage = this->tempDamage; +} + +void +LatticeFrameConcretePlasticStatus::printOutputAt(FILE *file, TimeStep *tStep) const +{ + LatticeMaterialStatus::printOutputAt(file, tStep); + + fprintf(file, "plasticStrains "); + for ( double s : this->plasticLatticeStrain ) { + fprintf(file, "% .8e ", s); + } + //fprintf(file, "kappad %.8e ", this->kappaD); + // fprintf(file, "damage %.8e ", this->damage); + fprintf(file, ", kappaD %.8e, damage %.8e \n", this->kappaD, this->damage); + +} + +void +LatticeFrameConcretePlasticStatus::saveContext(DataStream &stream, ContextMode mode) +// +// saves full information stored in this Status +// no temp variables stored +// +{ + LatticeMaterialStatus::saveContext(stream, mode); + + + + if ( !stream.write(& kappaD, 1) ) { + THROW_CIOERR(CIO_IOERR); + } + + if ( !stream.write(& damage, 1) ) { + THROW_CIOERR(CIO_IOERR); + } +} + + + +void +LatticeFrameConcretePlasticStatus::restoreContext(DataStream &stream, ContextMode mode) +// +// restores full information stored in stream to this Status +// +{ + LatticeMaterialStatus::restoreContext(stream, mode); + + if ( !stream.read(& kappaD, 1) ) { + THROW_CIOERR(CIO_IOERR); + } + + if ( !stream.read(& damage, 1) ) { + THROW_CIOERR(CIO_IOERR); + } +} +} diff --git a/src/sm/Materials/LatticeMaterials/latticeframeconcreteplastic.h b/src/sm/Materials/LatticeMaterials/latticeframeconcreteplastic.h new file mode 100644 index 000000000..bb75de887 --- /dev/null +++ b/src/sm/Materials/LatticeMaterials/latticeframeconcreteplastic.h @@ -0,0 +1,269 @@ +/* + * + * ##### ##### ###### ###### ### ### + * ## ## ## ## ## ## ## ### ## + * ## ## ## ## #### #### ## # ## + * ## ## ## ## ## ## ## ## + * ## ## ## ## ## ## ## ## + * ##### ##### ## ###### ## ## + * + * + * OOFEM : Object Oriented Finite Element Code + * + * Copyright (C) 1993 - 2022 Borek Patzak + * + * + * + * Czech Technical University, Faculty of Civil Engineering, + * Department of Structural Mechanics, 166 29 Prague, Czech Republic + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + */ + +#ifndef latticeframeconcreteplastic_h +#define latticeframeconcreteplastic_h + + +#include "latticestructuralmaterial.h" +#include "cltypes.h" +#include "randommaterialext.h" +#include "strainvector.h" +#include "stressvector.h" +#include "latticematstatus.h" +///@name Input fields for LatticeFrameConcretePlastic +//@{ +#define _IFT_LatticeFrameConcretePlastic_Name "LatticeFrameConcretePlastic" +#define _IFT_LatticeFrameConcretePlastic_talpha "talpha" +#define _IFT_LatticeFrameConcretePlastic_e "e" +#define _IFT_LatticeFrameConcretePlastic_n "n" +#define _IFT_LatticeFrameConcretePlastic_nx0 "nx0" +#define _IFT_LatticeFrameConcretePlastic_vy0 "vy0" +#define _IFT_LatticeFrameConcretePlastic_vz0 "vz0" +#define _IFT_LatticeFrameConcretePlastic_mx0 "mx0" +#define _IFT_LatticeFrameConcretePlastic_my0 "my0" +#define _IFT_LatticeFrameConcretePlastic_mz0 "mz0" +#define _IFT_LatticeFrameConcretePlastic_nx01 "nx01" +#define _IFT_LatticeFrameConcretePlastic_vy01 "vy01" +#define _IFT_LatticeFrameConcretePlastic_vz01 "vz01" +#define _IFT_LatticeFrameConcretePlastic_mx01 "mx01" +#define _IFT_LatticeFrameConcretePlastic_my01 "my01" +#define _IFT_LatticeFrameConcretePlastic_mz01 "mz01" +#define _IFT_LatticeFrameConcretePlastic_tol "tol" +#define _IFT_LatticeFrameConcretePlastic_iter "iter" +#define _IFT_LatticeFrameConcretePlastic_sub "sub" +#define _IFT_LatticeFrameConcretePlastic_plastic "plastic" +#define _IFT_LatticeFrameConcretePlastic_wu "wu" +#define _IFT_LatticeFrameConcretePlastic_wf "wf" +#define _IFT_LatticeFrameConcretePlastic_wftwo "wftwo" +#define _IFT_LatticeFrameConcretePlastic_qzero "qzero" +//@} + +namespace oofem { +/** + * This class implements associated Material Status to LatticeFrameConcretePlastic. + * @authors: Gumaa Abdelrhim, Peter Grassl + */ + +class LatticeFrameConcretePlasticStatus : public LatticeMaterialStatus +{ +public: + + enum state_flag_values { + LatticeFrameConcretePlastic_Elastic, + LatticeFrameConcretePlastic_Unloading, + LatticeFrameConcretePlastic_Plastic, + }; + enum LatticeFrameConcretePlastic_ReturnResult { + RR_NotConverged, + RR_Converged, + }; + +protected: + + double kappaD = 0.; + double tempKappaD = 0.; + double damage = 0.; + double tempDamage = 0.; + + + int tempReturnResult = LatticeFrameConcretePlasticStatus::RR_NotConverged; + +public: + + /// Constructor + LatticeFrameConcretePlasticStatus(int n, Domain *d, GaussPoint *g); + + + void printOutputAt(FILE *file, TimeStep *tStep) const override; + + const char *giveClassName() const override { return "LatticeFrameConcretePlasticStatus"; } + + void letTempReturnResultBe(const int result) { tempReturnResult = result; } + + int giveTempReturnResult() const { return tempReturnResult; } + + double giveKappaD() const { return this->kappaD; } + + double giveTempKappaD() const { return this->tempKappaD; } + + void setTempKappaD(double newKappa) { this->tempKappaD = newKappa; } + + double giveDamage() const { return this->damage; } + + double giveTempDamage() const { return this->tempDamage; } + + void setTempDamage(double newDamage) { this->tempDamage = newDamage; } + + void initTempStatus() override; + + void updateYourself(TimeStep *) override; + + void saveContext(DataStream &stream, ContextMode mode) override; + + void restoreContext(DataStream &stream, ContextMode mode) override; +}; + + +/** + * This class implements a concrete plasticity model for concrete for frame elements. + */ +class LatticeFrameConcretePlastic : public LatticeStructuralMaterial + +{ +protected: + + + + ///Normal modulus + double e; + + ///Ratio of shear and normal modulus + double nu; + + ///maximum axial force in x-axis x-axis nx0 + double nx0; + + double vy0; + + double vz0; + + ///maximum bending moment about x-axis mx0 + double mx0; + + ///maximum bending moment about x-axis my0 + double my0; + + ///maximum bending moment about x-axis mz0 + double mz0; + + ///maximum axial force in x-axis x-axis nx01 + double nx01; + + double vy01; + + double vz01; + + ///maximum bending moment about x-axis mx01 + double mx01; + + ///maximum bending moment about x-axis my01 + double my01; + + ///maximum bending moment about x-axis mz01 + double mz01; + + /// yield tolerance + double yieldTol; + + /// maximum number of iterations for stress return + double newtonIter; + + ///number Of SubIncrements + double numberOfSubIncrements; + + ///plastic flag + double plasticFlag; + + ///wu + double wu; + + ///wf + double wf; + + ///wftwo + double wftwo; + + ///qzero + double qzero; + + + enum LatticeFrameConcretePlastic_ReturnResult { RR_NotConverged, RR_Converged, RR_Unknown, RR_known }; + + double initialYieldStress = 0.; + +public: + LatticeFrameConcretePlastic(int n, Domain *d) : LatticeStructuralMaterial(n, d) { }; + + FloatArrayF< 6 >computeFVector(const FloatArrayF< 6 > &stress, const FloatArrayF< 6 > &k, GaussPoint *gp, TimeStep *tStep) const; + + FloatMatrixF< 6, 6 >computeDMMatrix(const FloatArrayF< 6 > &stress, const FloatArrayF< 6 > &k, GaussPoint *gp, TimeStep *tStep) const; + + FloatArrayF< 6 >giveThermalDilatationVector(GaussPoint *gp, TimeStep *tStep) const override; + + FloatArrayF< 6 >giveLatticeStrain(GaussPoint *gp, TimeStep *tStep) const; + + int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *atTime) override; + virtual FloatArrayF< 6 >giveStrain(GaussPoint *gp, TimeStep *tStep) const; + + FloatArrayF< 6 >performPlasticityReturn(GaussPoint *gp, const FloatArrayF< 6 > &strain, TimeStep *tStep) const; + + void performRegularReturn(FloatArrayF< 6 > &stress, LatticeFrameConcretePlastic_ReturnResult, const FloatArrayF< 6 > &k, double yieldValue, GaussPoint *gp, TimeStep *tStep) const; + + void giveSwitches(IntArray &answer, int k); + + double computeYieldValue(const FloatArrayF< 6 > &stress, const FloatArrayF< 6 > &k, GaussPoint *gp, TimeStep *tStep) const; + + FloatMatrixF< 7, 7 >computeJacobian(const FloatArrayF< 6 > &stress, const FloatArrayF< 6 > &k, const double deltaLambda, GaussPoint *gp, TimeStep *tStep) const; + + IntArray checkStatus(const FloatArrayF< 6 > &stress, IntArray &k, GaussPoint *gp, TimeStep *tStep) const; + + FloatArrayF< 6 >giveFrameForces3d(const FloatArrayF< 6 > &strain, GaussPoint *gp, TimeStep *tStep) override; + + const char *giveInputRecordName() const override { return _IFT_LatticeFrameConcretePlastic_Name; } + + const char *giveClassName() const override { return "LatticeFrameConcretePlastic"; } + + void initializeFrom(InputRecord &ir) override; + + bool isCharacteristicMtrxSymmetric(MatResponseMode rMode) const override { return false; } + + Interface *giveInterface(InterfaceType) override; + + FloatMatrixF< 6, 6 >give3dFrameStiffnessMatrix(MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep) const override; + + bool hasMaterialModeCapability(MaterialMode mode) const override; + + MaterialStatus *CreateStatus(GaussPoint *gp) const override; + + MaterialStatus *giveStatus(GaussPoint *gp) const override; + +protected: + const FloatArray checkStatus(FloatArrayF< 6 >f, GaussPoint *pPoint, TimeStep *pStep) const; + IntArray checkStatus(const FloatArrayF< 6 > &stress, IntArray &k) const; + const FloatArrayF< 6 >checkTransition(const FloatArrayF< 6 > &stress, GaussPoint *gp, TimeStep *tStep) const; +}; +} // end namespace oofem + +#endif diff --git a/src/sm/Materials/LatticeMaterials/latticeframesteelplastic.C b/src/sm/Materials/LatticeMaterials/latticeframesteelplastic.C index 62a74616b..02f6315a4 100644 --- a/src/sm/Materials/LatticeMaterials/latticeframesteelplastic.C +++ b/src/sm/Materials/LatticeMaterials/latticeframesteelplastic.C @@ -206,7 +206,7 @@ LatticeFrameSteelPlastic::computeDMMatrix(const FloatArrayF< 4 > &stress, GaussP dm.at(4, 1) = 0; dm.at(4, 2) = 0; dm.at(4, 3) = 0; - dm.at(4, 3) = 2. / pow(this->mz0, 2.); + dm.at(4, 4) = 2. / pow(this->mz0, 2.); return dm; } diff --git a/src/sm/Materials/LatticeMaterials/latticematstatus.C b/src/sm/Materials/LatticeMaterials/latticematstatus.C index e5f3e1c8f..55a734652 100644 --- a/src/sm/Materials/LatticeMaterials/latticematstatus.C +++ b/src/sm/Materials/LatticeMaterials/latticematstatus.C @@ -54,6 +54,8 @@ LatticeMaterialStatus :: initTempStatus() this->tempLatticeStrain = this->latticeStrain; + this->tempInternalForces = this->internalForces; + this->tempLatticeStress = this->latticeStress; this->tempReducedLatticeStrain = this->reducedLatticeStrain; @@ -78,6 +80,8 @@ LatticeMaterialStatus :: updateYourself(TimeStep *atTime) { MaterialStatus :: updateYourself(atTime); + this->internalForces = this->tempInternalForces; + this->latticeStress = this->tempLatticeStress; this->latticeStrain = this->tempLatticeStrain; @@ -146,6 +150,10 @@ LatticeMaterialStatus :: saveContext(DataStream &stream, ContextMode mode) THROW_CIOERR(iores); } + if ( ( iores = internalForces.storeYourself(stream) ) != CIO_OK ) { + THROW_CIOERR(iores); + } + if ( ( iores = latticeStrain.storeYourself(stream) ) != CIO_OK ) { THROW_CIOERR(iores); } @@ -162,6 +170,10 @@ LatticeMaterialStatus :: saveContext(DataStream &stream, ContextMode mode) THROW_CIOERR(iores); } + if ( ( iores = internalForces.storeYourself(stream) ) != CIO_OK ) { + THROW_CIOERR(iores); + } + if ( !stream.write(le) ) { THROW_CIOERR(CIO_IOERR); } @@ -194,6 +206,10 @@ LatticeMaterialStatus :: restoreContext(DataStream &stream, ContextMode mode) THROW_CIOERR(iores); } + if ( ( iores = internalForces.restoreYourself(stream) ) != CIO_OK ) { + THROW_CIOERR(iores); + } + if ( ( iores = latticeStress.restoreYourself(stream) ) != CIO_OK ) { THROW_CIOERR(iores); } @@ -210,6 +226,11 @@ LatticeMaterialStatus :: restoreContext(DataStream &stream, ContextMode mode) THROW_CIOERR(iores); } + if ( ( iores = internalForces.restoreYourself(stream) ) != CIO_OK ) { + THROW_CIOERR(iores); + } + + if ( !stream.read(le) ) { THROW_CIOERR(CIO_IOERR); } @@ -225,5 +246,7 @@ LatticeMaterialStatus :: restoreContext(DataStream &stream, ContextMode mode) if ( !stream.read(crackFlag) ) { THROW_CIOERR(CIO_IOERR); } + + } } // end namespace oofem diff --git a/src/sm/Materials/LatticeMaterials/latticematstatus.h b/src/sm/Materials/LatticeMaterials/latticematstatus.h index 125136184..b5093fc4b 100644 --- a/src/sm/Materials/LatticeMaterials/latticematstatus.h +++ b/src/sm/Materials/LatticeMaterials/latticematstatus.h @@ -78,6 +78,10 @@ class LatticeMaterialStatus : public MaterialStatus, public RandomMaterialStatus /// Non-equilibrated plastic lattice strain FloatArrayF< 6 >oldPlasticLatticeStrain; + /// Equilibrated Internal Forces + FloatArrayF< 12 >internalForces; + ///Non-equilibrated Internal Forces + FloatArrayF< 12 >tempInternalForces; /// Equilibriated damage lattice strain FloatArrayF< 6 >damageLatticeStrain; @@ -101,6 +105,12 @@ class LatticeMaterialStatus : public MaterialStatus, public RandomMaterialStatus ///Increment of dissipation double deltaDissipation = 0.; + // ///Internal lForces + // double internalForces=0; + + // ///temp Internal Forces + // double tempInternalForces=0; + /// Non-equilibrated increment of dissipation double tempDeltaDissipation = 0.; @@ -139,8 +149,10 @@ class LatticeMaterialStatus : public MaterialStatus, public RandomMaterialStatus /// Returns lattice strain. const FloatArrayF< 6 > &giveLatticeStrain() const { return this->latticeStrain; } + + /// Returns lattice strain. - const FloatArrayF< 6 > &giveTempLatticeStrain() const { return this->tempLatticeStress; } + const FloatArrayF< 6 > &giveTempLatticeStrain() const { return this->tempLatticeStrain; } /// Returns reduced lattice strain. const FloatArrayF< 6 > &giveReducedLatticeStrain() const { return reducedLatticeStrain; } @@ -162,9 +174,16 @@ class LatticeMaterialStatus : public MaterialStatus, public RandomMaterialStatus /// Returns temp lattice stress. const FloatArrayF< 6 > &giveTempLatticeStress() const { return this->tempLatticeStress; } + /// Returns temp lattice stress. + const FloatArrayF< 12 > &giveTempInternalForces() const { return this->tempInternalForces; } + /// Returns temp damage lattice strain. const FloatArrayF< 6 > &giveTempDamageLatticeStrain() const { return this->tempDamageLatticeStrain; } + /// Returns temp damage lattice strain. + const FloatArrayF< 12 > &giveInternalForces() const { return this->internalForces; } + + /// Assigns the temp value of lattice strain. void letTempLatticeStrainBe(const FloatArrayF< 6 > &v) { this->tempLatticeStrain = v; } @@ -177,6 +196,9 @@ class LatticeMaterialStatus : public MaterialStatus, public RandomMaterialStatus /// Assigns the temp value of lattice stress. void letTempLatticeStressBe(const FloatArrayF< 6 > &v) { this->tempLatticeStress = v; } + /// Assigns the temp value of lattice stress. + void letTempInternalForcesBe(const FloatArrayF< 12 > &v) { this->tempInternalForces = v; } + /// Assigns the temp value of damage lattice strain. void letTempDamageLatticeStrainBe(const FloatArrayF< 6 > &v) { this->tempDamageLatticeStrain = v; } diff --git a/tests/sm/latticeframe3dnl.in b/tests/sm/latticeframe3dnl.in new file mode 100644 index 000000000..085905c77 --- /dev/null +++ b/tests/sm/latticeframe3dnl.in @@ -0,0 +1,42 @@ +latticeframe3dnl.out +Buckling3D. Model +NonlinearStatic nmsteps 1 nsteps 1 nmodules 1 +nsteps 100 stiffmode 0 rtolv 1.e-3 maxiter 20000 initialguess 1 controlmode 1 +errorcheck +#vtkxml tstep_step 1 primvars 1 1 cellvars 8 138 139 140 47 46 45 141 52 domain_all element_all +domain 3dLattice +OutputManager +ndofman 9 nelem 8 ncrosssect 1 nmat 1 nbc 2 nic 0 nltf 2 nset 0 +node 1 coords 3 0.0 0.0 0.0 bc 6 1 1 1 1 1 1 +node 2 coords 3 1.1775 0.0 0.0 +node 3 coords 3 2.355 0.0 0.0 +node 4 coords 3 3.5325 0.0 0.0 +node 5 coords 3 4.71 0.0 0.0 +node 6 coords 3 5.8875 0.0 0.0 +node 7 coords 3 7.065 0.0 0.0 +node 8 coords 3 8.2425 0.0 0.0 +node 9 coords 3 9.42 0.0 0.0 bc 6 0 0 0 0 0 2 +# +latticeframe3dnl 1 nodes 2 1 2 zaxis 3 0 0 1 crossSect 1 mat 1 +latticeframe3dnl 2 nodes 2 2 3 zaxis 3 0 0 1 crossSect 1 mat 1 +latticeframe3dnl 3 nodes 2 3 4 zaxis 3 0 0 1 crossSect 1 mat 1 +latticeframe3dnl 4 nodes 2 4 5 zaxis 3 0 0 1 crossSect 1 mat 1 +latticeframe3dnl 5 nodes 2 5 6 zaxis 3 0 0 1 crossSect 1 mat 1 +latticeframe3dnl 6 nodes 2 6 7 zaxis 3 0 0 1 crossSect 1 mat 1 +latticeframe3dnl 7 nodes 2 7 8 zaxis 3 0 0 1 crossSect 1 mat 1 +latticeframe3dnl 8 nodes 2 8 9 zaxis 3 0 0 1 crossSect 1 mat 1 +# +latticecs 1 material 1 iy 7.857e-9 iz 7.857e-9 ik 1.6e-8 shearCoeff 0.83 area 3.1416e-4 +latticeframeelastic 1 d 1. E 2.1e11 n 0.3 talpha 0 +BoundaryCondition 1 loadTimeFunction 1 prescribedvalue 0.0 +BoundaryCondition 2 loadTimeFunction 2 prescribedvalue 1. +ConstantFunction 1 f(t) 1 +PiecewiseLinFunction 2 t 3 -1. 0. 100 f(t) 3 0. 0. .1 +#%BEGIN_CHECK% tolerance 1.e-6 +#NODE tStep 100 number 1 dof 1 unknown d value 0.0 +#NODE tStep 100 number 1 dof 2 unknown d value 0.0 +#NODE tStep 100 number 5 dof 1 unknown d value -1.89407969e-03 +#NODE tStep 100 number 5 dof 2 unknown d value 1.16549434e-01 +#NODE tStep 100 number 9 dof 1 unknown d value -1.53218429e-02 +#NODE tStep 100 number 9 dof 2 unknown d value 4.65912214e-01 +#%END_CHECK% diff --git a/tests/sm/latticeframesteelplastic1.in b/tests/sm/latticeframesteelplastic1.in index 10a2a8c5e..ad2600f99 100644 --- a/tests/sm/latticeframesteelplastic1.in +++ b/tests/sm/latticeframesteelplastic1.in @@ -1,9 +1,10 @@ latticeframesteelplastic1.out Cantilever test with latticeframesteelplastic material model -NonLinearStatic nmsteps 1 nsteps 20 contextOutputStep 100 nmodules 0 +NonLinearStatic nmsteps 1 nsteps 20 contextOutputStep 100 nmodules 1 nsteps 20 rtolv 1.e-0 maxiter 2000 stiffMode 1 controllmode 1 refloadmode 1 minsteplength 1.e-3 ddm 2 3 6 ddv 1 2 ddltf 2 +errorcheck domain 3dLattice -OutputManager tstep_all dofman_all element_all +OutputManager tstep_all ndofman 3 nelem 2 ncrosssect 1 nmat 1 nbc 1 nic 0 nltf 2 node 1 coords 3 0.000000e+00 0.000000e+00 0.000000e+00 bc 6 1 1 1 1 1 1 node 2 coords 3 0.500000e+00 0.000000e+00 0.000000e+00 @@ -15,20 +16,7 @@ latticeframesteelplastic 1 d 1. E 1. n 0.15 talpha 0 nx0 1 mx0 1 my0 1 mz0 1 tol BoundaryCondition 1 loadTimeFunction 1 prescribedvalue 0.0 ConstantFunction 1 f(t) 1. PiecewiseLinFunction 2 nPoints 2 t 2 0 20. f(t) 2 0. 2 -#PiecewiseLinFunction 3 nPoints 2 t 2 0 10. f(t) 2 0. 1 -#PiecewiseLinFunction 2 nPoints 5 t 5 0. 60 61 62 300. f(t) 5 0. 40. 18.8 40. 100 #%BEGIN_CHECK% -#REACTION number 1 dof 1 -#REACTION number 1 dof 2 -#REACTION number 1 dof 3 -#REACTION number 1 dof 4 -#REACTION number 1 dof 5 -#REACTION number 1 dof 6 -#DOFMAN number 3 dof 1 type d -#DOFMAN number 3 dof 2 type d -#DOFMAN number 3 dof 3 type d -#DOFMAN number 3 dof 4 type d -#DOFMAN number 3 dof 5 type d -#DOFMAN number 3 dof 6 type d -#TIME +#NODE tStep 20 number 3 dof 2 unknown d value 1.90000118e+01 tolerance 1.e-4 +#NODE tStep 20 number 3 dof 6 unknown d value 3.80000237e+01 tolerance 1.e-4 #%END_CHECK% diff --git a/tests/sm/latticeframesteelplastic2.in b/tests/sm/latticeframesteelplastic2.in index 763ed2d4c..4e82f2353 100644 --- a/tests/sm/latticeframesteelplastic2.in +++ b/tests/sm/latticeframesteelplastic2.in @@ -1,9 +1,11 @@ latticeframesteelplastic2.out Toi 1993 plastic frame -NonLinearStatic nmsteps 1 nsteps 1 contextOutputStep 100 nmodules 0 -nsteps 10 rtolv 1.e-6 maxiter 100 stiffMode 1 controllmode 1 refloadmode 1 minsteplength 1.e-3 ddm 2 5 1 ddv 1 2e-1 ddltf 2 +NonLinearStatic nmsteps 1 nsteps 1 contextOutputStep 100 nmodules 1 +nsteps 10 rtolv 1.e-6 maxiter 100 stiffMode 1 controllmode 1 refloadmode 1 minsteplength 1.e-3 ddm 2 5 1 ddv 1 2e-1 ddltf 2 +#vtkxml tstep_step 1 primvars 1 1 domain_all element_all +errorcheck domain 3dLattice -OutputManager tstep_all dofman_all element_all +OutputManager ndofman 13 nelem 12 ncrosssect 1 nmat 1 nbc 1 nic 0 nltf 2 node 1 coords 3 0.000000e+00 0.000000e+00 0.000000e+00 bc 6 1 1 1 1 1 1 node 2 coords 3 0.000000e+00 250.000000e+00 0.000000e+00 @@ -36,17 +38,14 @@ BoundaryCondition 1 loadTimeFunction 1 prescribedvalue 0.0 ConstantFunction 1 f(t) 1. PiecewiseLinFunction 2 nPoints 2 t 2 0 10. f(t) 2 0. 10 #%BEGIN_CHECK% -#REACTION number 1 dof 1 -#REACTION number 1 dof 2 -#REACTION number 1 dof 3 -#REACTION number 1 dof 4 -#REACTION number 1 dof 5 -#REACTION number 1 dof 6 -#DOFMAN number 5 dof 1 type d -#DOFMAN number 5 dof 2 type d -#DOFMAN number 5 dof 3 type d -#DOFMAN number 5 dof 4 type d -#DOFMAN number 5 dof 5 type d -#DOFMAN number 5 dof 6 type d +#REACTION tStep 10 number 1 dof 1 value -7.13485006e+04 tolerance 1.e-3 +#REACTION tStep 10 number 1 dof 2 value -6.49983931e+04 tolerance 1.e-3 +#REACTION tStep 10 number 1 dof 6 value 3.87420101e+07 tolerance 1. +#REACTION tStep 10 number 13 dof 1 value -7.11033739e+04 tolerance 1.e-3 +#REACTION tStep 10 number 13 dof 2 value 6.49983991e+04 tolerance 1.e-3 +#REACTION tStep 10 number 13 dof 6 value 3.87113692e+07 tolerance 1. +#NODE tStep 10 number 5 dof 1 unknown d value 9.00000895e+00 tolerance 1.e-5 +#NODE tStep 10 number 5 dof 2 unknown d value 1.36093742e-01 tolerance 1.e-5 +#NODE tStep 10 number 5 dof 6 unknown d value -5.26397500e-03 tolerance 1.e-6 #TIME #%END_CHECK%