fem » ele

coupled.h

00001 /*************************************************************************************
00002  * MechSys - A C++ library to simulate (Continuum) Mechanical Systems                *
00003  * Copyright (C) 2005 Dorival de Moraes Pedroso <dorival.pedroso at gmail.com>       *
00004  * Copyright (C) 2005 Raul Dario Durand Farfan  <raul.durand at gmail.com>           *
00005  *                                                                                   *
00006  * This file is part of MechSys.                                                     *
00007  *                                                                                   *
00008  * MechSys is free software; you can redistribute it and/or modify it under the      *
00009  * terms of the GNU General Public License as published by the Free Software         *
00010  * Foundation; either version 2 of the License, or (at your option) any later        *
00011  * version.                                                                          *
00012  *                                                                                   *
00013  * MechSys is distributed in the hope that it will be useful, but WITHOUT ANY        *
00014  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A   *
00015  * PARTICULAR PURPOSE. See the GNU General Public License for more details.          *
00016  *                                                                                   *
00017  * You should have received a copy of the GNU General Public License along with      *
00018  * MechSys; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, *
00019  * Fifth Floor, Boston, MA 02110-1301, USA                                           *
00020  *************************************************************************************/
00021 
00022 #ifndef FEM_COUPLED_H
00023 #define FEM_COUPLED_H
00024 
00025 #ifndef REAL 
00026     #define REAL double
00027 #endif
00028 
00029 #include <sstream>
00030 
00031 #include "models/coupled/coupledmodel.h"
00032 #include "fem/ele/element.h"
00033 #include "fem/ele/equilibelem.h"
00034 #include "fem/ele/flowelem.h"
00035 #include "linalg/laexpr.h"
00036 #include "tensors/tensors.h"
00037 #include "tensors/functions.h"
00038 #include "util/numstreams.h"
00039 
00040 using Tensors::Stress_p_q_S_sin3th; // Mean and deviatoric Cambridge stress invariants
00041 using Tensors::Strain_Ev_Ed;        // Volumetric and deviatoric strain invariants
00042 namespace FEM
00043 {
00044     
00045 
00046 class Coupled: public virtual EquilibElem, public virtual FlowElem
00047 {
00048     static const int NDIM=3;
00049     static const int NSTRESSCOMPS=6;
00050 public:
00051     // Constructor
00052     Coupled();
00053     // Destructor
00054     virtual ~Coupled() {};
00055     // Methods
00056     void Deactivate();
00057     void ReAllocateModel(String const & ModelName, Array<REAL> const & ModelPrms, Array<Array<REAL> > const & AIniData);
00058     bool IsEssential(String const & DOFName) const;
00059     void CalcFaceNodalValues(String            const & FaceDOFName  , // In:  Face DOF Name, ex.: "tx, ty, tz"
00060                              REAL              const   FaceDOFValue , // In:  Face DOF Value, ex.: 10 kPa, 20 kPa
00061                              Array<FEM::Node*> const & APtrFaceNodes, // In:  Array of ptrs to face nodes
00062                              String                  & NodalDOFName , // Out: Nodal DOF Name, ex.: "fx, fy, fz"
00063                              LinAlg::Vector<REAL>    & NodalValues    // Out: Vector of nodal values
00064     ) const;
00065     void Stiffness       (LinAlg::Matrix<REAL> & Ke, Array<size_t> & EqMap) const;
00066     void CouplingMatrix1 (LinAlg::Matrix<REAL> & L1, Array<size_t> & EqMapEquilib, Array<size_t> & EqMapFlow) const; 
00067     void CouplingMatrix2 (LinAlg::Matrix<REAL> & L2, Array<size_t> & EqMapEquilib, Array<size_t> & EqMapFlow) const; 
00068     void MassMatrix      (LinAlg::Matrix<REAL> & M,  Array<size_t> & EqMap) const;
00069     void Permeability    (LinAlg::Matrix<REAL> & H,  Array<size_t> & EqMap) const; 
00070     // Create a structure to hold the DOFs information of Node [j_node] inside Element [i_ele]
00071     // These DOFs are the only ones that Element [i_ele] uses => NOT the global ones
00072     void NodalDOFs(int iNode, Array<FEM::Node::DOFVarsStruct*> & DOFs) const;
00073     void BackupState();
00074     void UpdateState(REAL TimeInc); // Read nodal values _IncEssential and write to _IncNatural
00075     void RestoreState();
00076     void SetProperties(Array<REAL> const & EleProps) { _unit_weight=EleProps[0]; }
00077     String OutCenter(bool PrintCaptionOnly) const;
00078     void OutNodes(LinAlg::Matrix<REAL> & Values, Array<String> & Labels) const;
00079     // Functions to assemble DAS matrices
00080     size_t nOrder0Matrices() const { return 1; };
00081     size_t nOrder1Matrices() const { return 4; };
00082     void      RHSVector(size_t index, LinAlg::Vector<REAL> & V, Array<size_t> & Map) const;
00083     void   Order0Matrix(size_t index, LinAlg::Matrix<REAL> & V, Array<size_t> & RowsMap, Array<size_t> & ColsMap) const;
00084     void   Order1Matrix(size_t index, LinAlg::Matrix<REAL> & M, Array<size_t> & RowsMap, Array<size_t> & ColsMap) const;
00085     // Output
00086     REAL OutScalar1()             const; // Pore-pressure
00087     REAL OutScalar2()             const; // Ed
00088     void OutTensor1(String & Str) const; // Stress
00089     void OutTensor2(String & Str) const; // Strain
00090 private:
00091     // Data
00092     Array<CoupledModel*> _a_model;
00093     REAL                 _unit_weight;
00094     // Auxiliar Data
00095     mutable Array<Tensors::Tensor2> _a_d_ep; // Elastoplastic vectors (stress-suction) obtained from Stiffness stored temporarily to be used later by CouplingMatrix1
00096     // Methods
00097     void _set_node_vars(int iNode);
00098     void _calc_Nu_matrix(LinAlg::Vector<REAL> & Shape, LinAlg::Matrix<REAL> & Nu) const;
00099     void _calc_initial_internal_forces();
00100     void _calc_initial_internal_pore_pressures();
00101 }; // class Coupled
00102 
00103 // ========================================================================== Implementation
00104 
00105 // Constructor
00106 Coupled::Coupled() // {{{
00107 {
00108     _a_d_ep.resize(_n_int_pts); //  Resizing array of d_ep vectors 
00109 } // }}}
00110 
00111 inline void Coupled::_set_node_vars(int iNode) // {{{
00112 {
00113     // Add Degree of Freedom to a node (Essential, Natural) 
00114     _connects[iNode]->AddDOF(DUX,DFX);
00115     _connects[iNode]->AddDOF(DUY,DFY);
00116     _connects[iNode]->AddDOF(DUZ,DFZ);
00117     _connects[iNode]->AddDOF(DWP,DWD);
00118 } // }}}
00119 
00120 void Coupled::Deactivate() // {{{
00121 {
00122     if (_is_active==false) return; // the element was already inactive
00123     _is_active = false; 
00124 
00125     //dereferencing nodes:
00126     for (int i_node=0; i_node<_n_nodes; i_node++)
00127     {
00128         _connects[i_node]->Refs()--;
00129     }
00130 
00131     //verifying if element is in the boundary of excavation
00132     bool in_boundary = false;
00133     for (int i_node=0; i_node<_n_nodes; i_node++)
00134         if (_connects[i_node]->Refs()>0) 
00135         {
00136             in_boundary=true;
00137             break;
00138         }
00139 
00140     //calculate internal forces for element on boudary
00141     if (in_boundary)
00142     {
00143         LinAlg::Vector<REAL> F(NDIM*_n_nodes);
00144         F.SetValues(0.0);
00145 
00146         // Allocate entities used for every integration point
00147         LinAlg::Matrix<REAL> derivs; // size = NumLocalCoords(ex.: r,s,t) x _n_nodes
00148         LinAlg::Matrix<REAL> J;      // Jacobian matrix
00149         LinAlg::Matrix<REAL> B;      // strain-displacement matrix
00150         LinAlg::Vector<REAL> Sig;    // stress tensor
00151         
00152         for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
00153         {
00154             // Temporary Integration Points
00155             REAL r = _a_int_pts[i_ip].r;
00156             REAL s = _a_int_pts[i_ip].s;
00157             REAL t = _a_int_pts[i_ip].t;
00158             REAL w = _a_int_pts[i_ip].w;
00159 
00160             Derivs(r,s,t, derivs);        // Calculate Derivatives of Shape functions w.r.t local coordinate system
00161             Jacobian(derivs, J);          // Calculate J (Jacobian) matrix for i_ip Integration Point
00162             B_Matrix(derivs,J, B);        // Calculate B matrix for i_ip Integration Point
00163 
00164             // Get tensor for accumulate stress
00165             Tensors::Tensor2 const & tsrSig = _a_model[i_ip]->Sig();
00166             Copy(tsrSig, Sig);
00167 
00168             // Calculate internal force vector;
00169             F += trn(B)*Sig*det(J)*w;
00170         }   
00171 
00172         // adding to boundary values of node if Refs()>0
00173         for (int i_node=0; i_node<_n_nodes; i_node++)
00174             if (_connects[i_node]->Refs()>0) 
00175             {
00176                 _connects[i_node]->DOFVar(DFX).NaturalBry += F(i_node*NDIM  );
00177                 _connects[i_node]->DOFVar(DFY).NaturalBry += F(i_node*NDIM+1);
00178                 _connects[i_node]->DOFVar(DFZ).NaturalBry += F(i_node*NDIM+2);
00179                 
00180                 _connects[i_node]->DOFVar(DWP).IsEssenPresc = true;
00181                 // REAL pwp = _connects[i_node]->DOFVar(DWP).EssentialVal; 
00182                 _connects[i_node]->DOFVar(DWP).EssentialBry = -_connects[i_node]->DOFVar(DWP).EssentialVal;
00183                 _connects[i_node]->DOFVar(DWP).EssenPrescPersist = true; // that solution is ugly and other solution must be formulated
00184             }
00185             else // eliminate dofs from nodes with no references
00186             {
00187                 _connects[i_node]->DOFVars().resize(0);
00188             }
00189 
00190         //delete IP information 
00191     }
00192 } 
00193 
00194 void Coupled::ReAllocateModel(String const & ModelName, Array<REAL> const & ModelPrms, Array<Array<REAL> > const & AIniData) // {{{
00195 {
00196     // Check size of AIniData
00197     if (!(AIniData.size()==1 || static_cast<int>(AIniData.size())==_n_int_pts))
00198         throw new Fatal(_("Coupled::ReAllocateModel: Array of array of initial data must have size==1 or size== %d \n"), _n_int_pts);
00199 
00200     // If pointers to model was not already defined => No model was allocated
00201     if (_a_model.size()==0)
00202     {
00203         // Resize the array of model pointers
00204         _a_model.resize(_n_int_pts);
00205 
00206         // Loop along integration points
00207         for (int i=0; i<_n_int_pts; ++i)
00208         {
00209             // Allocate a new model and set parameters
00210             if (AIniData.size()==1) _a_model[i] = AllocCoupledModel(ModelName, ModelPrms, AIniData[0]);
00211             else                    _a_model[i] = AllocCoupledModel(ModelName, ModelPrms, AIniData[i]);
00212         }
00213 
00214         // Calculate initial internal forces
00215         _calc_initial_internal_forces();
00216         // Calculate initial internal pore-pressures
00217         _calc_initial_internal_pore_pressures();
00218     }
00219     else // Model objects already defined (allocated)
00220     {
00221         // Loop along integration points
00222         for (int i=0; i<_n_int_pts; ++i)
00223         {
00224             if (_a_model[i]->Name()!=ModelName) // Do re-allocate
00225             {
00226                 // Delete old model
00227                 delete _a_model[i];
00228 
00229                 // Allocate a new model and set parameters
00230                 if (AIniData.size()==1) _a_model[i] = AllocCoupledModel(ModelName, ModelPrms, AIniData[0]);
00231                 else                    _a_model[i] = AllocCoupledModel(ModelName, ModelPrms, AIniData[i]);
00232             }
00233             // else: Do NOT re-allocate => Model not changed (don't do anything)
00234         }
00235     }
00236 } // }}}
00237 
00238 inline bool Coupled::IsEssential(String const & DOFName) const // {{{
00239 {
00240     if (DOFName==DUX || DOFName==DUY || DOFName==DUZ || DOFName==DWP) return true;
00241     else return false;
00242 } // }}}
00243 
00244 inline void Coupled::CalcFaceNodalValues(String            const & FaceDOFName  ,         // In:  Face DOF Name, ex.: "tx, ty, tz" {{{
00245                                          REAL              const   FaceDOFValue ,         // In:  Face DOF Value, ex.: 10 kPa, 20 kPa
00246                                          Array<FEM::Node*> const & APtrFaceNodes,         // In:  Array of ptrs to face nodes
00247                                          String                  & NodalDOFName ,         // Out: Nodal DOF Name, ex.: "fx, fy, fz"
00248                                          LinAlg::Vector<REAL>    & NodalValues  ) const   // Out: Vector of nodal values
00249 {
00250     if (FaceDOFName==DTX || FaceDOFName==DTY || FaceDOFName==DTZ || FaceDOFName==DFWD)
00251     {
00252         if (FaceDOFName==DTX) NodalDOFName=DFX;
00253         if (FaceDOFName==DTY) NodalDOFName=DFY;
00254         if (FaceDOFName==DTZ) NodalDOFName=DFZ;
00255         if (FaceDOFName==DFWD) NodalDOFName=DWD;
00256         _distrib_val_to_face_nodal_vals(APtrFaceNodes, FaceDOFValue, NodalValues);
00257     }
00258     else
00259     {
00260         std::ostringstream oss; oss << "Face nodes coordinates:\n";
00261         for (size_t i_node=0; i_node<APtrFaceNodes.size(); ++i_node)
00262             oss << "X=" << APtrFaceNodes[i_node]->X() << ", Y=" << APtrFaceNodes[i_node]->Y() << ", Z=" << APtrFaceNodes[i_node]->Z() << std::endl;
00263         throw new Fatal(_("Coupled::CalcFaceNodalValues: This method must only be called for FaceDOFName< %s > equal to tx, ty, tz or tq\n %s \n"), FaceDOFName.c_str(), oss.str().c_str());
00264     }
00265 } // }}}
00266 
00267 inline void Coupled::Stiffness(LinAlg::Matrix<REAL> & Ke, Array<size_t> & EqMap) const // {{{
00268 {
00269     // Resize Ke
00270     Ke.Resize(NDIM*_n_nodes, NDIM*_n_nodes); // sum(Bt*D*B*det(J)*w)
00271     Ke.SetValues(0.0);
00272 
00273     // Allocate entities used for every integration point
00274     LinAlg::Matrix<REAL> derivs; // size = NumLocalCoords(ex.: r,s,t) x _n_nodes
00275     LinAlg::Matrix<REAL> J;      // Jacobian matrix
00276     LinAlg::Matrix<REAL> B;      // strain-displacement matrix
00277     Tensors::Tensor4     tsrD;   // Fourth order constitutive tensor
00278     LinAlg::Matrix<REAL> D;      // Constitutive matrix
00279 
00280     // Loop along integration points
00281     for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
00282     {
00283         // Temporary Integration Points
00284         REAL r = _a_int_pts[i_ip].r;
00285         REAL s = _a_int_pts[i_ip].s;
00286         REAL t = _a_int_pts[i_ip].t;
00287         REAL w = _a_int_pts[i_ip].w;
00288 
00289         Derivs(r,s,t, derivs);        // Calculate Derivatives of Shape functions w.r.t local coordinate system
00290         Jacobian(derivs, J);          // Calculate J (Jacobian) matrix for i_ip Integration Point
00291         B_Matrix(derivs,J, B);        // Calculate B matrix for i_ip Integration Point
00292 
00293         // Constitutive tensor & elastoplastic vector d_ep to be used later by CouplingMatrix1
00294         _a_model[i_ip]->TgStiffness(tsrD, _a_d_ep[i_ip]); Copy(tsrD, D);
00295 
00296         // Calculate Tangent Stiffness
00297         Ke += trn(B)*D*B*det(J)*w;
00298     }
00299 
00300     //Mounting a map of positions from Ke to Global
00301     int idx_Ke=0; // position (idx) inside Ke matrix
00302     EqMap.resize(Ke.Rows()); // size=Ke.Rows()=Ke.Cols()
00303     for (int i_node=0; i_node<_n_nodes; ++i_node)
00304     {
00305         // Fill map of Ke position to K position of DOFs components
00306         EqMap[idx_Ke++] = _connects[i_node]->DOFVar(DUX).EqID; 
00307         EqMap[idx_Ke++] = _connects[i_node]->DOFVar(DUY).EqID; 
00308         EqMap[idx_Ke++] = _connects[i_node]->DOFVar(DUZ).EqID; 
00309     }
00310 } // }}} 
00311 
00312 inline void Coupled::CouplingMatrix1(LinAlg::Matrix<REAL> & L1, Array<size_t> & EqMapEquilib, Array<size_t> & EqMapFlow) const // {{{
00313 {
00314     //  
00315     //   Coupling Matrix L1:  ( n_dim*_n_nodes x _n_nodes )
00316     //   ==================================================
00317     //       
00318     //                  /    T      T                            /    T        T            
00319     //        [L1]   =  | [B]  * {m} * phi(Sr) * {N}  * dV   -   | [B]  * {dep} * {N}  * dV  // coupled unsaturated
00320     //                  /    u                      p            /    u             p      
00321     //       
00322     //                  /    T      T             
00323     //        [L1]   =  | [B]  * {m} * {N}  * dV                                             // coupled saturated
00324     //                  /    u            p  
00325     //           
00326     //            T                   
00327     //         {m}   = [ 1 1 1 0 0 0 ]
00328     
00329     // Resize L1
00330     L1.Resize(NDIM*_n_nodes, 1*_n_nodes); 
00331     L1.SetValues(0.0);
00332 
00333     // Allocate entities used for every integration point
00334     LinAlg::Matrix<REAL> derivs; // size = NumLocalCoords(ex.: r,s,t) x _n_nodes
00335     LinAlg::Vector<REAL> shape;  // size = _n_nodes
00336     LinAlg::Matrix<REAL> J;      // Jacobian matrix
00337     LinAlg::Matrix<REAL> B;      // strain-displacement matrix
00338     LinAlg::Vector<REAL> d_ep;   // elasto-plastic constitutive vector 
00339     LinAlg::Vector<REAL> m(6);
00340     m = 1, 1, 1, 0, 0, 0;
00341 
00342     // Loop along integration points
00343     for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
00344     {
00345         // Temporary Integration Points
00346         REAL r = _a_int_pts[i_ip].r;
00347         REAL s = _a_int_pts[i_ip].s;
00348         REAL t = _a_int_pts[i_ip].t;
00349         REAL w = _a_int_pts[i_ip].w;
00350 
00351         Shape(r,s,t, shape);                // Calculate Np vector (equal to shape functions vector)
00352         Derivs(r,s,t, derivs);              // Calculate Derivatives of Shape functions w.r.t local coordinate system
00353         Jacobian(derivs, J);                // Calculate J (Jacobian) matrix for i_ip Integration Point
00354         B_Matrix(derivs,J, B);              // Calculate B matrix for i_ip Integration Point
00355         REAL det_J = det(J);
00356         REAL phi   = _a_model[i_ip]->phi(); // Calculate phi(Sr) for i_ip Integration Point
00357 
00358         Copy(_a_d_ep[i_ip], d_ep);
00359         L1 += trn(B)*m*trn(shape)*phi*det_J*w - trn(B)*d_ep*trn(shape)*det_J*w;
00360     }
00361     
00362     //Mounting a map of positions from L1 to Global
00363     int idx_L1=0; // position (idx) inside L1 matrix
00364     EqMapEquilib.resize(L1.Rows()); // size=L1.Rows()
00365     for (int i_node=0; i_node<_n_nodes; ++i_node)
00366     {
00367         // Fill map of L1 position to Global position of DOFs components
00368         EqMapEquilib[idx_L1++] = _connects[i_node]->DOFVar(DUX).EqID; 
00369         EqMapEquilib[idx_L1++] = _connects[i_node]->DOFVar(DUY).EqID; 
00370         EqMapEquilib[idx_L1++] = _connects[i_node]->DOFVar(DUZ).EqID; 
00371     }
00372     
00373     idx_L1=0; // position (idx) inside L1 matrix
00374     EqMapFlow.resize(L1.Cols()); // size=L1.Cols()
00375     for (int i_node=0; i_node<_n_nodes; ++i_node)
00376     {
00377         // Fill map of L1 position to Global position of DOFs components
00378         EqMapFlow[idx_L1++] = _connects[i_node]->DOFVar(DWP).EqID; 
00379     }
00380     
00381 } // }}}
00382 
00383 inline void Coupled::CouplingMatrix2(LinAlg::Matrix<REAL> & L2, Array<size_t> & EqMapFlow, Array<size_t> & EqMapEquilib) const // {{{
00384 {
00385     //  
00386     //   Coupling Matrix L2:  ( _n_nodes x n_dim*_n_nodes )
00387     //   =================================================
00388     //       
00389     //                  /                T      
00390     //         [L2]  =  | {N}  * Sr * {m} * [B]  * dV     // coupled unsaturated
00391     //                  /    p                u  
00392     //
00393     //   p: pore-pressure
00394     //   u: displacement
00395     //
00396     //      T
00397     //   {m} = [ 1 1 1 0 0 0 ]
00398     //           
00399     
00400     // Resize L2
00401     L2.Resize(1*_n_nodes, NDIM*_n_nodes); // sum(Bt*m*Np*det(J)*w)
00402     L2.SetValues(0.0);
00403 
00404     // Allocate entities used for every integration point
00405     LinAlg::Matrix<REAL> derivs; // size = NumLocalCoords(ex.: r,s,t) x _n_nodes
00406     LinAlg::Vector<REAL> shape;  // size = _n_nodes
00407     LinAlg::Matrix<REAL> J;      // Jacobian matrix
00408     LinAlg::Matrix<REAL> B;      // strain-displacement matrix
00409     LinAlg::Vector<REAL> m(6);
00410     m = 1, 1, 1, 0, 0, 0;
00411 
00412     // Loop along integration points
00413     for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
00414     {
00415         // Temporary Integration Points
00416         REAL r = _a_int_pts[i_ip].r;
00417         REAL s = _a_int_pts[i_ip].s;
00418         REAL t = _a_int_pts[i_ip].t;
00419         REAL w = _a_int_pts[i_ip].w;
00420 
00421         Shape(r,s,t, shape);              // Calculate Np vector (equal to shape functions vector)
00422         Derivs(r,s,t, derivs);            // Calculate Derivatives of Shape functions w.r.t local coordinate system
00423         Jacobian(derivs, J);              // Calculate J (Jacobian) matrix for i_ip Integration Point
00424         B_Matrix(derivs,J, B);            // Calculate B matrix for i_ip Integration Point
00425         REAL Sr = _a_model[i_ip]->Sr();   // Degree of saturation
00426 
00427         L2 += shape*Sr*trn(m)*B*det(J)*w;
00428     }
00429     
00430     //Mounting a map of positions from L2 to Global
00431     
00432     int idx_L2=0; // position (idx) inside L2 matrix
00433     EqMapFlow.resize(L2.Rows()); // size=L2.Rows()
00434     for (int i_node=0; i_node<_n_nodes; ++i_node)
00435     {
00436         // Fill map of L2 position to Global position of DOFs components
00437         EqMapFlow[idx_L2++] = _connects[i_node]->DOFVar(DWP).EqID; 
00438     }
00439 
00440     idx_L2=0; // position (idx) inside L2 matrix
00441     EqMapEquilib.resize(L2.Cols()); // size=L2.Cols()
00442     for (int i_node=0; i_node<_n_nodes; ++i_node)
00443     {
00444         // Fill map of L2 position to Global position of DOFs components
00445         EqMapEquilib[idx_L2++] = _connects[i_node]->DOFVar(DUX).EqID; 
00446         EqMapEquilib[idx_L2++] = _connects[i_node]->DOFVar(DUY).EqID; 
00447         EqMapEquilib[idx_L2++] = _connects[i_node]->DOFVar(DUZ).EqID; 
00448     }
00449 } // }}}
00450 
00451 inline void Coupled::MassMatrix(LinAlg::Matrix<REAL> & M, Array<size_t> & EqMap) const // {{{
00452 {
00453     //  
00454     //   Mass Matrix M: ( _n_nodes x _n_nodes );
00455     //   ============================
00456     //       
00457     //                 /                        T 
00458     //         [M]  =  | {N}  * n * dSr_dp * {N}  * dV     // coupled unsaturated
00459     //                 /    p                  p  
00460     //
00461     //   p: pore-pressure
00462     //           
00463     
00464     // Resize M
00465     M.Resize(1*_n_nodes,1*_n_nodes); 
00466     M.SetValues(0.0);
00467 
00468     // Allocate entities used for every integration point
00469     LinAlg::Matrix<REAL> derivs; // size = NumLocalCoords(ex.: r,s,t) x _n_nodes
00470     LinAlg::Vector<REAL> shape;  // size = _n_nodes
00471     LinAlg::Matrix<REAL> J;      // Jacobian matrix
00472     LinAlg::Vector<REAL> m(6);
00473     m = 1, 1, 1, 0, 0, 0;
00474 
00475     // Loop along integration points
00476     for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
00477     {
00478         // Temporary Integration Points
00479         REAL r = _a_int_pts[i_ip].r;
00480         REAL s = _a_int_pts[i_ip].s;
00481         REAL t = _a_int_pts[i_ip].t;
00482         REAL w = _a_int_pts[i_ip].w;
00483 
00484         Shape(r,s,t, shape);              // Calculate Np vector (equal to shape functions vector)
00485         Derivs(r,s,t, derivs);            // Calculate Derivatives of Shape functions w.r.t local coordinate system
00486         Jacobian(derivs, J);              // Calculate J (Jacobian) matrix for i_ip Integration Point
00487         REAL n_times_dSr_dPp = _a_model[i_ip]->n_times_dSr_dPp();
00488 
00489         // Calculate the M matrix
00490         M += -shape*trn(shape)*n_times_dSr_dPp*det(J)*w;
00491     }
00492     
00493     //Mounting a map of positions from M to Global
00494     int idx_S=0; // position (idx) inside M matrix
00495     EqMap.resize(M.Rows()); // size=M.Rows==M.Cols()
00496     for (int i_node=0; i_node<_n_nodes; ++i_node)
00497     {
00498         // Fill map of M position to Global position of DOFs components
00499         EqMap[idx_S++] = _connects[i_node]->DOFVar(DWP).EqID; 
00500     }
00501     
00502 } // }}}
00503 
00504 inline void Coupled::Permeability(LinAlg::Matrix<REAL> & He, Array<size_t> & EqMap) const // {{{
00505 {
00506     //  
00507     //   Permeability Matrix K:
00508     //   ============================
00509     //       
00510     //                  /    T                   
00511     //         [P]   =  | [B]  * [K] * [B]  * dV
00512     //                  /    p            p  
00513     //
00514     //   OBS: [K] = [k]/gamaW        
00515     
00516     // Resize He
00517     He.Resize(_n_nodes, _n_nodes); // sum(Bpt*k*Bp*det(J)*w)
00518     He.SetValues(0.0);
00519 
00520     // Allocate entities used for every integration point
00521     LinAlg::Matrix<REAL> derivs; // size = NumLocalCoords(ex.: r,s,t) x _n_nodes
00522     LinAlg::Vector<REAL> shape;  // size = _n_nodes
00523     LinAlg::Matrix<REAL> J;      // Jacobian matrix
00524     LinAlg::Matrix<REAL> Bp;     // pore-pressure - gradient matrix
00525     Tensor2              tsrK;   // permeability tensor
00526     LinAlg::Matrix<REAL> K;      // permeability matrix
00527 
00528     // Loop along integration points
00529     for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
00530     {
00531         // Temporary Integration Points
00532         REAL r = _a_int_pts[i_ip].r;
00533         REAL s = _a_int_pts[i_ip].s;
00534         REAL t = _a_int_pts[i_ip].t;
00535         REAL w = _a_int_pts[i_ip].w;
00536         Shape(r,s,t, shape);              // Calculate Np vector (equal to shape functions vector)
00537         Derivs(r,s,t, derivs);            // Calculate Derivatives of Shape functions w.r.t local coordinate system
00538         Jacobian(derivs, J);              // Calculate J (Jacobian) matrix for i_ip Integration Point
00539         Bp_Matrix(derivs,J, Bp);            // Calculate Bp matrix for i_ip Integration Point
00540 
00541         // Calculate Tangent Permeability
00542         _a_model[i_ip]->TgPermeability(tsrK); Copy(tsrK, K);
00543         REAL gammaW = _a_model[i_ip]->GammaW();
00544 
00545         // Calculate the permeability matrix
00546         He += -trn(Bp)*K*Bp*det(J)*w/gammaW;
00547     }
00548     
00549     //Mounting a map of positions from He to Global
00550     int idx_He=0; // position (idx) inside He matrix
00551     EqMap.resize(He.Rows()); // size=He.Rows()=He.Cols()
00552     for (int i_node=0; i_node<_n_nodes; ++i_node)
00553     {
00554         // Fill map of He position to Global position of DOFs components
00555         EqMap[idx_He++] = _connects[i_node]->DOFVar(DWP).EqID; 
00556     }
00557 } // }}}
00558 
00559 inline void Coupled::_calc_Nu_matrix(LinAlg::Vector<REAL> & Shape, //{{{
00560                                     LinAlg::Matrix<REAL> & Nu//Nu must be (3 x 3*_n_nodes);
00561                                     ) const
00562 {
00563     Nu.Resize(NDIM, NDIM*_n_nodes);
00564     for(int i=0; i<_n_nodes; i++)
00565     {
00566         Nu(0,0 + 3*i) = Shape(i);  Nu(0,1 + 3*i) = 0;         Nu(0,2 + 3*i) = 0;    
00567         Nu(1,0 + 3*i) = 0;         Nu(1,1 + 3*i) = Shape(i);  Nu(1,2 + 3*i) = 0;    
00568         Nu(2,0 + 3*i) = 0;         Nu(2,1 + 3*i) = 0;         Nu(2,2 + 3*i) = Shape(i); 
00569     }
00570 } //}}}
00571 
00572 inline void Coupled::NodalDOFs(int iNode, Array<FEM::Node::DOFVarsStruct*> & DOFs) const // {{{
00573 {
00574     // Resize DOFs array: 3 displacement + 1 pore-pressure
00575     DOFs.resize(4);
00576 
00577     // Set DOFs array with pointers to DOFVarsStruct inside this element Node [iNode]
00578     DOFs[0] = &_connects[iNode]->DOFVar(DUX);
00579     DOFs[1] = &_connects[iNode]->DOFVar(DUY);
00580     DOFs[2] = &_connects[iNode]->DOFVar(DUZ);
00581     DOFs[3] = &_connects[iNode]->DOFVar(DWP);
00582 
00583 } // }}}
00584 
00585 inline void Coupled::BackupState() // {{{
00586 
00587 {
00588     for (int i=0; i<_n_int_pts; ++i)
00589         _a_model[i]->BackupState();
00590 } // }}}
00591 
00592 inline void Coupled::UpdateState(REAL TimeInc) // Read nodal values _IncEssential and write to _IncNatural  {{{
00593 {
00594     //   Saturated:
00595     //   ==========
00596     //
00597     //   Internal force integration:                                   Internal volume integration:                                                                 
00598     //                                                                     
00599     //            int    /    T                                                 int     /    T                 /    T
00600     //         {dF}   =  | [B]  * [dSig] * dV                                {dQ}   =  -| [N]  * dn * dV   +   | [B]  * h * {v} * dV
00601     //                   /                                                              /    p                 /    p               
00602     //                                                                         
00603     //         (dSig is related to total stress increment)                   h = TimeInc                                            
00604     //           
00605     //         {dSig} = {dSig'}  +  dp * {m}  
00606     //
00607     //                 T         
00608     //         dp = [N]  * {dP} 
00609     //  
00610     //   Unsaturated:
00611     //   ============
00612     //
00613     //   Internal force integration:                                   Internal volume integration:                                         
00614     //                                                                                                                                  
00615     //            int    /    T                                           
00616     //         {dF}   =  | [B]  * {dSig} * dV                                   int     /    T                     /    T
00617     //                   /                                                   {dQ}   = - | [N]  * d(nSr) * dV   +   | [B]  * h * {v} * dV
00618     //                                                                                  /    p                     /    p               
00619     //         (dSig is related to total stress increment)                     
00620     //                                                                       h = TimeInc                                                
00621     //         {dSig} = {dSig'}  +  (phi_before*dp + p*dphi)*{m};
00622     //  
00623     //                 T                        T       
00624     //         dp = [N]  * {dP}         p = [N]  * {P}
00625     //
00626 
00627     // Allocate (local/element) vectors
00628     LinAlg::Vector<REAL> dU(NDIM*_n_nodes); // Delta disp. of  this element
00629     LinAlg::Vector<REAL> dP(_n_nodes);      // Delta pore-pres sure of this element
00630     LinAlg::Vector<REAL>  P(_n_nodes);       // Total pore-pressure of this element
00631     
00632     // Loop along nodes
00633     for (int i_node=0; i_node<_n_nodes; ++i_node)
00634     {
00635         // Assemble (local/element) displacements vector. OBS.: NDIM=3
00636         dU(i_node*NDIM  ) = _connects[i_node]->DOFVar(DUX)._IncEssenVal;
00637         dU(i_node*NDIM+1) = _connects[i_node]->DOFVar(DUY)._IncEssenVal;
00638         dU(i_node*NDIM+2) = _connects[i_node]->DOFVar(DUZ)._IncEssenVal;
00639         dP(i_node)        = _connects[i_node]->DOFVar(DWP)._IncEssenVal;
00640         P(i_node)         = _connects[i_node]->DOFVar(DWP).EssentialVal;
00641     }
00642      
00643     // Allocate (local/element) internal force vector
00644     LinAlg::Vector<REAL> dF(NDIM*_n_nodes); // Delta internal force of this element
00645     dF.SetValues(0.0);
00646     
00647     // Allocate (local/element) internal volume vector
00648     LinAlg::Vector<REAL> dQ(1*_n_nodes); // Delta internal volumes of this element
00649     dQ.SetValues(0.0);
00650 
00651     // Allocate entities used for every integration point
00652     LinAlg::Matrix<REAL> derivs;  // size = NumLocalCoords(ex.: r,s,t) x _n_nodes
00653     LinAlg::Vector<REAL> shape;   // size = _n_nodes
00654     LinAlg::Matrix<REAL> J;       // Jacobian matrix 
00655     LinAlg::Matrix<REAL> B;       // strain-displacedment matrix
00656     LinAlg::Matrix<REAL> Bp;      // pore-pressure - gradient matrix
00657     Tensors::Tensor2     tsrDEps; // Strain tensor  tsrDEps(NSTRESSCOMPS=6) = B(NSTRESSCOMPS=6,NDIM*_n_nodes) * dU(NDIM*_n_nodes)
00658     Tensors::Tensor2     tsrDSig; // Stress tensor
00659     Tensors::Tensor1     tsrVel;  // tensor for velocity
00660     Tensors::Tensor1     tsrGrad; // tensor for gradient
00661     LinAlg::Vector<REAL> DEps;    // Strain vector in Mandel's notation 
00662     LinAlg::Vector<REAL> DSig;    // Stress vector in Mandel's notation 
00663     LinAlg::Vector<REAL> Vel;     // vector for velocity
00664     LinAlg::Vector<REAL> Grad;    // vector for gradient
00665     LinAlg::Vector<REAL> m(6);
00666     m = 1, 1, 1, 0, 0, 0;
00667 
00668     // Loop along integration points
00669     for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
00670     {
00671         // Temporary Integration Points
00672         REAL r = _a_int_pts[i_ip].r;
00673         REAL s = _a_int_pts[i_ip].s;
00674         REAL t = _a_int_pts[i_ip].t;
00675         REAL w = _a_int_pts[i_ip].w;
00676 
00677         Shape(r,s,t, shape);              // Calculate Np vector (equal to shape functions vector)
00678         Derivs(r,s,t, derivs);            // Calculate Derivatives of Shape functions w.r.t local coordinate system
00679         Jacobian(derivs, J);              // Calculate J (Jacobian) matrix for i_ip Integration Point
00680         B_Matrix(derivs, J, B);           // Calculate B  matrix for i_ip Integration Point
00681         Bp_Matrix(derivs, J, Bp);         // Calculate Bp matrix for i_ip Integration Point
00682         REAL det_J = J.Det();
00683 
00684         // Calculate internal forces ==================================================
00685 
00686         // Calculate a tensor for the increments of strain
00687         DEps = B*dU; Copy(DEps, tsrDEps); // tsrDEps = B*dU
00688 
00689         // Update model
00690         REAL          p = _a_model[i_ip]->Pp();   // total value of pore-pressure
00691         REAL         dp = trn(shape)*dP;          // increment of pore-pressure
00692         REAL phi_before = _a_model[i_ip]->phi();
00693         REAL DnSr; // will be updated by StressUpdate method
00694         _a_model[i_ip]->StressUpdate(tsrDEps, dp, tsrDSig, DnSr); Copy(tsrDSig, DSig);
00695         _a_model[i_ip]->FlowUpdate(dp);
00696         // Adding increments of pore-pressure to DSig for positive values of total pore-pressure
00697         // For        Sig' = Sig   - phi*uw*m   in which: Sig = total streses, Sig' = (net/constitutive) stresses
00698         // then:      Sig  = Sig'  + phi*uw*m             uw  = p (pore-pressure)
00699         //           dSig  = dSig' + (phi*dp + p*dphi)*m
00700         REAL phi_after = _a_model[i_ip]->phi();
00701         REAL dphi = phi_after - phi_before;
00702         DSig += (phi_before*dp + p*dphi)*m;
00703 
00704         // Calculate internal force vector
00705         dF += trn(B)*DSig*det_J*w;
00706 
00707         // Calculate internal volumes =================================================
00708 
00709         REAL gammaW = _a_model[i_ip]->GammaW();
00710         Grad = Bp*P/gammaW; Copy(Grad, tsrGrad);
00711         _a_model[i_ip]->FlowVelocity(tsrGrad, tsrVel); Copy(tsrVel, Vel);
00712 
00713         // Calculate internal volumes
00714         dQ += -shape*DnSr*det_J*w + trn(Bp)*TimeInc*Vel*det_J*w;
00715     }
00716 
00717     // Update nodal _IncNaturVals
00718     for (int i_node=0; i_node<_n_nodes; ++i_node)
00719     {
00720         // Assemble (local/element) displacements vector. OBS.: NDIM=3
00721         _connects[i_node]->DOFVar(DFX)._IncNaturVal += dF(i_node*NDIM  );
00722         _connects[i_node]->DOFVar(DFY)._IncNaturVal += dF(i_node*NDIM+1);
00723         _connects[i_node]->DOFVar(DFZ)._IncNaturVal += dF(i_node*NDIM+2);
00724         _connects[i_node]->DOFVar(DWD) ._IncNaturVal += dQ(i_node);
00725 
00726         // Update natural vars state
00727         _connects[i_node]->DOFVar(DFX).NaturalVal += dF(i_node*NDIM  );
00728         _connects[i_node]->DOFVar(DFY).NaturalVal += dF(i_node*NDIM+1);
00729         _connects[i_node]->DOFVar(DFZ).NaturalVal += dF(i_node*NDIM+2);
00730         _connects[i_node]->DOFVar(DWD) .NaturalVal += dQ(i_node);
00731     }
00732 
00733 } // }}}
00734 
00735 inline void Coupled::RestoreState() // {{{
00736 {
00737     for (int i=0; i<_n_int_pts; ++i)
00738         _a_model[i]->RestoreState();
00739 } // }}}
00740 
00741 inline String Coupled::OutCenter(bool PrintCaptionOnly=false) const // {{{
00742 {
00743     // Auxiliar variables
00744     std::ostringstream oss;
00745 
00746     // Number of state values
00747     int n_int_state_vals = _a_model[0]->nInternalStateValues();
00748     
00749     // Print caption
00750     if (PrintCaptionOnly)
00751     {
00752         // Stress and strains
00753         oss << _8s() <<  "p" << _8s() << "q"  << _8s() << "th" << _8s() << "Ev"  << _8s() << "Ed";
00754         oss << _8s() << "Sx" << _8s() << "Sy" << _8s() << "Sz" << _8s() << "Sxy" << _8s() << "Syz" << _8s() << "Sxz";
00755         oss << _8s() << "Ex" << _8s() << "Ey" << _8s() << "Ez" << _8s() << "Exy" << _8s() << "Eyz" << _8s() << "Exz";
00756         oss << _8s() << "Pp";
00757 
00758         // Internal state names
00759         Array<String> str_state_names;   _a_model[0]->InternalStateNames(str_state_names);
00760         for (int i=0; i<n_int_state_vals; ++i)
00761             oss << _8s() << str_state_names[i];
00762         oss << std::endl;
00763     }
00764     else
00765     {
00766         // Stress, strains and internal state values evaluated at the center of the element
00767         Tensors::Tensor2 sig_cen(0.0);
00768         Tensors::Tensor2 eps_cen(0.0);
00769         REAL             ppr_cen=0.0; // pore-pressure
00770         Array<REAL>      int_state_vals_cen;   int_state_vals_cen.assign(n_int_state_vals,0.0);
00771 
00772         // Loop over integration points
00773         for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
00774         {
00775             // Stress and strains
00776             sig_cen += _a_model[i_ip]->Sig();
00777             eps_cen += _a_model[i_ip]->Eps();
00778             ppr_cen += _a_model[i_ip]->Pp();
00779 
00780             // Internal state values
00781             Array<REAL> int_state_vals;    _a_model[i_ip]->InternalStateValues(int_state_vals);
00782             for (int j=0; j<n_int_state_vals; ++j)
00783                 int_state_vals_cen[j] += int_state_vals[j];
00784         }
00785         
00786         // Average stress and strains
00787         sig_cen = sig_cen / _n_int_pts;
00788         eps_cen = eps_cen / _n_int_pts;
00789         ppr_cen = ppr_cen / _n_int_pts;
00790         
00791         // Average internal state values
00792         for (int j=0; j<n_int_state_vals; ++j)
00793             int_state_vals_cen[j] = int_state_vals_cen[j] / _n_int_pts;;
00794 
00795         // Calculate stress invariants
00796         REAL              p,q,sin3th;
00797         Tensors::Tensor2  S;
00798         Stress_p_q_S_sin3th(sig_cen,p,q,S,sin3th);
00799 
00800         // Calculate strain invariants
00801         REAL     Ev,Ed;
00802         Strain_Ev_Ed(eps_cen,Ev,Ed);
00803 
00804         // Output
00805         REAL sq2 = sqrt(2.0);
00806         oss << _8s()<< p          << _8s()<< q          << _8s()<< sin3th << _8s()<< Ev*100.0       << _8s()<< Ed*100.0;
00807         oss << _8s()<< sig_cen(0) << _8s()<< sig_cen(1) << _8s()<< sig_cen(2)         << _8s()<< sig_cen(3)/sq2 << _8s()<< sig_cen(4)/sq2 << _8s()<< sig_cen(5)/sq2;
00808         oss << _8s()<< eps_cen(0) << _8s()<< eps_cen(1) << _8s()<< eps_cen(2)         << _8s()<< eps_cen(3)/sq2 << _8s()<< eps_cen(4)/sq2 << _8s()<< eps_cen(5)/sq2;
00809         oss << _8s()<< ppr_cen;
00810 
00811         // Internal state values
00812         for (int j=0; j<n_int_state_vals; ++j)
00813             oss << _8s()<< int_state_vals_cen[j];
00814         oss << std::endl;
00815 
00816     }
00817 
00818     return oss.str(); 
00819 } // }}}
00820 
00821 inline void Coupled::OutNodes(LinAlg::Matrix<REAL> & Values, Array<String> & Labels) const // {{{
00822 {
00823     int const DATA_COMPS=20;
00824     Values.Resize(_n_nodes,DATA_COMPS);
00825     Labels.resize(DATA_COMPS);
00826     Labels[ 0] = DUX;  Labels[ 1] = DUY;  Labels[ 2] = DUZ;  Labels[ 3] = DWP;
00827     Labels[ 4] = DFX;  Labels[ 5] = DFY;  Labels[ 6] = DFZ;  Labels[ 7] = DWD;
00828     Labels[ 8] = "Ex"; Labels[ 9] = "Ey"; Labels[10] = "Ez"; Labels[11] = "Exy"; Labels[12] = "Eyz"; Labels[13] = "Exz";
00829     Labels[14] = "Sx"; Labels[15] = "Sy"; Labels[16] = "Sz"; Labels[17] = "Sxy"; Labels[18] = "Syz"; Labels[19] = "Sxz";
00830     for (int i_node=0; i_node<_n_nodes; i_node++)
00831     {
00832         Values(i_node,0) = _connects[i_node]->DOFVar(DUX).EssentialVal;
00833         Values(i_node,1) = _connects[i_node]->DOFVar(DUY).EssentialVal;
00834         Values(i_node,2) = _connects[i_node]->DOFVar(DUZ).EssentialVal;
00835         Values(i_node,3) = _connects[i_node]->DOFVar(DWP).EssentialVal;
00836         Values(i_node,4) = _connects[i_node]->DOFVar(DFX).NaturalVal;
00837         Values(i_node,5) = _connects[i_node]->DOFVar(DFY).NaturalVal;
00838         Values(i_node,6) = _connects[i_node]->DOFVar(DFZ).NaturalVal;
00839         Values(i_node,7) = _connects[i_node]->DOFVar(DWD).NaturalVal;
00840     }
00841 
00842     //Extrapolation
00843     LinAlg::Vector<REAL> ip_values(_n_int_pts);
00844     LinAlg::Vector<REAL> nodal_values(_n_nodes);
00845     
00846     // Strains
00847     for (int i_comp=0; i_comp<NSTRESSCOMPS; i_comp++)
00848     {
00849         for (int j_ip=0; j_ip<_n_int_pts; j_ip++)
00850         {
00851             Tensors::Tensor2 const & eps = _a_model[j_ip]->Eps();
00852             ip_values(j_ip) = eps(i_comp); //getting IP values
00853         }
00854         _extrapolate(ip_values, nodal_values);
00855         for (int j_node=0; j_node<_n_nodes; j_node++)
00856             Values(j_node,i_comp+8) = nodal_values(j_node);
00857     }
00858     
00859     // Stresses
00860     for (int i_comp=0; i_comp<NSTRESSCOMPS; i_comp++)
00861     {
00862         for (int j_ip=0; j_ip<_n_int_pts; j_ip++)
00863         {
00864             Tensors::Tensor2 const & sig = _a_model[j_ip]->Sig();
00865             ip_values(j_ip) = sig(i_comp); //getting IP values
00866         }
00867         _extrapolate(ip_values, nodal_values);
00868         for (int j_node=0; j_node<_n_nodes; j_node++)
00869             Values(j_node,i_comp+14) = nodal_values(j_node);
00870     }
00871 } // }}}
00872 
00873 inline void Coupled::_calc_initial_internal_forces() // {{{
00874 {
00875     // Allocate (local/element) internal force vector
00876     LinAlg::Vector<REAL> F(NDIM*_n_nodes);
00877     F.SetValues(0.0);
00878 
00879     // Allocate entities used for every integration point
00880     LinAlg::Matrix<REAL> derivs;  // size = NumLocalCoords(ex.: r,s,t) x _n_nodes
00881     LinAlg::Matrix<REAL> J;       // Jacobian matrix
00882     LinAlg::Matrix<REAL> B;       // strain-displacement matrix
00883     Tensors::Tensor2     tsrSig;  // Stress tensor 
00884     LinAlg::Vector<REAL> Sig;     // Stress vector in Mandel's notation 
00885 
00886     // Loop along integration points
00887     for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
00888     {
00889         // Temporary Integration Points
00890         REAL r = _a_int_pts[i_ip].r;
00891         REAL s = _a_int_pts[i_ip].s;
00892         REAL t = _a_int_pts[i_ip].t;
00893         REAL w = _a_int_pts[i_ip].w;
00894 
00895         Derivs(r,s,t, derivs);        // Calculate Derivatives of Shape functions w.r.t local coordinate system
00896         Jacobian(derivs, J);          // Calculate J (Jacobian) matrix for i_ip Integration Point
00897         B_Matrix(derivs, J, B);       // Calculate B matrix for i_ip Integration Point
00898 
00899         REAL phi = _a_model[i_ip]->phi();
00900         tsrSig   = _a_model[i_ip]->Sig() + phi*_a_model[i_ip]->Pp()*Tensors::I;
00901         Copy(tsrSig, Sig);
00902 
00903         // Calculate internal force vector;
00904         F += trn(B)*Sig*det(J)*w;
00905     }
00906 
00907     // Update nodal NaturVals
00908     for (int i_node=0; i_node<_n_nodes; ++i_node)
00909     {
00910         // Assemble (local/element) displacements vector. OBS.: NDIM=3
00911         _connects[i_node]->DOFVar(DFX).NaturalVal += F(i_node*NDIM  ); // NaturalVal must be set to zero during AddDOF routine
00912         _connects[i_node]->DOFVar(DFY).NaturalVal += F(i_node*NDIM+1);
00913         _connects[i_node]->DOFVar(DFZ).NaturalVal += F(i_node*NDIM+2);
00914     }
00915 } // }}}
00916 
00917 inline void Coupled::_calc_initial_internal_pore_pressures() // {{{
00918 {
00919     // Here must be done an extrapolation form pore-pressures values into integration
00920     // points to nodal points
00921     // For simplicity, here is considered the mean value of integration points
00922     //
00923     // Here is considered an initial stationary state, it is the reason ...
00924     // ... to leave the calculus of initital internal volumes.
00925     
00926     // Allocate (local/element) internal force vector
00927     LinAlg::Vector<REAL> P(_n_nodes);
00928     P.SetValues(0.0);
00929 
00930     REAL mean_pp_value=0;
00931     // Loop along integration points
00932     for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
00933     {
00934         mean_pp_value+= _a_model[i_ip]->Pp()/_n_int_pts;
00935     }
00936 
00937     for (int i_node=0; i_node<_n_int_pts; ++i_node)
00938     {
00939         P(i_node) = mean_pp_value;
00940     }
00941 
00942     // Update nodal NaturVals
00943     for (int i_node=0; i_node<_n_nodes; ++i_node)
00944     {
00945         // Assemble (local/element) pore-pressures vector
00946         int node_refs = _connects[i_node]->Refs();
00947         _connects[i_node]->DOFVar(DWP).EssentialVal += P(i_node)/node_refs; // NaturalVal must be set to zero during AddDOF routine
00948         //_connects[i_node]->DOFVar("q").NaturalVal += Q(i_node); // NaturalVal must be set to zero during AddDOF routine
00949     }
00950 } // }}}
00951 
00952 inline void Coupled::RHSVector(size_t index, LinAlg::Vector<REAL> & V, Array<size_t> & Map) const // {{{
00953 {
00954     assert(index == 0);
00955     V.Resize(_n_nodes);
00956     
00957     // Loop along nodes
00958     for (int i_node=0; i_node<_n_nodes; ++i_node)
00959     {
00960         // Assemble (local/element) total pore-presure vector. 
00961         V(i_node) = _connects[i_node]->DOFVar(DWP).EssentialVal;
00962     }
00963     //Mounting a map of positions f
00964     int idx=0; // position (idx) inside V vector
00965     Map.resize(V.Size()); 
00966     for (int i_node=0; i_node<_n_nodes; ++i_node)
00967     {
00968         Map[idx++] = _connects[i_node]->DOFVar(DWP).EqID; 
00969     }
00970 } // }}}
00971 
00972 inline void Coupled::Order0Matrix(size_t index, LinAlg::Matrix<REAL> & M, Array<size_t> & RowsMap, Array<size_t> & ColsMap) const // {{{
00973 {
00974     assert(index == 0);
00975     Permeability(M, RowsMap);
00976     ColsMap.resize(RowsMap.size());
00977     for (size_t i=0; i<RowsMap.size(); i++) ColsMap[i] = RowsMap[i];
00978 } // }}}
00979 
00980 inline void Coupled::Order1Matrix(size_t index, LinAlg::Matrix<REAL> & M, Array<size_t> & RowsMap, Array<size_t> & ColsMap) const // {{{
00981 {
00982     assert(index < 4);
00983     switch(index)
00984     {
00985         case 0: 
00986             Stiffness(M, RowsMap);
00987             ColsMap.resize(RowsMap.size());
00988             for (size_t i=0; i<RowsMap.size(); i++) ColsMap[i] = RowsMap[i];
00989             break;
00990         case 1: 
00991             CouplingMatrix1(M, RowsMap, ColsMap);
00992             break;
00993         case 2: 
00994             CouplingMatrix2(M, RowsMap, ColsMap);
00995             break;
00996         case 3: 
00997             MassMatrix(M, RowsMap);
00998             ColsMap.resize(RowsMap.size());
00999             for (size_t i=0; i<RowsMap.size(); i++) ColsMap[i] = RowsMap[i];
01000             break;
01001     }
01002 } // }}}
01003 
01004 inline REAL Coupled::OutScalar1() const // {{{ Pore-pressure
01005 {
01006     // Pore-pressure evaluated at the center or the element
01007     REAL ppr=0.0;
01008 
01009     // Loop over integration points
01010     for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
01011         ppr += _a_model[i_ip]->Pp();
01012     
01013     // Average pore-pressure
01014     ppr = ppr / _n_int_pts;
01015 
01016     // Output
01017     return ppr;
01018 
01019 } // }}}
01020 
01021 inline REAL Coupled::OutScalar2() const // {{{ Ed
01022 {
01023     // Strains evaluated at the center of the element
01024     Tensors::Tensor2 e(0.0);
01025 
01026     // Loop over integration points
01027     for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
01028         e += 100.0*_a_model[i_ip]->Eps();
01029     
01030     // Average strains
01031     e = e / _n_int_pts;
01032 
01033     // Calculate strain invariants
01034     REAL   Ev,Ed;
01035     Strain_Ev_Ed(e,Ev,Ed);
01036 
01037     // Output
01038     return Ed;
01039 
01040 } // }}}
01041 
01042 inline void Coupled::OutTensor1(String & Str) const // {{{ Stress
01043 {
01044     // Stress evaluated at the center of the element
01045     Tensors::Tensor2 s(0.0);
01046 
01047     // Loop over integration points
01048     for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
01049         s += _a_model[i_ip]->Sig();
01050     
01051     // Average stress
01052     s = s / _n_int_pts;
01053 
01054     // Output
01055     REAL sq2 = sqrt(2.0);
01056     Str.Printf(_(" %e %e %e  %e %e %e  %e %e %e "), s(0),s(3)/sq2,s(5)/sq2,  s(3)/sq2,s(1),s(4)/sq2,  s(5)/sq2,s(4)/sq2,s(2));
01057 
01058 } // }}}
01059 
01060 inline void Coupled::OutTensor2(String & Str) const // {{{ Strain
01061 {
01062     // Strains evaluated at the center of the element
01063     Tensors::Tensor2 e(0.0);
01064 
01065     // Loop over integration points
01066     for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
01067         e += 100.0*_a_model[i_ip]->Eps();
01068     
01069     // Average strains
01070     e = e / _n_int_pts;
01071 
01072     // Output
01073     REAL sq2 = sqrt(2.0);
01074     Str.Printf(_(" %e %e %e  %e %e %e  %e %e %e "), e(0),e(3)/sq2,e(5)/sq2,  e(3)/sq2,e(1),e(4)/sq2,  e(5)/sq2,e(4)/sq2,e(2));
01075 
01076 } // }}}
01077 
01078 // Register the DOF information of CoupledElement into DOFInfoMap
01079 int CoupledDOFInfoRegister() // {{{
01080 {
01081     // Temporaries
01082     DOFInfo D; 
01083     DOFInfo D_flow; 
01084     D_flow = DOFInfoMap["FlowElem"];
01085 
01086     // EquilibElem DOF info
01087     D      = DOFInfoMap["EquilibElem"];
01088 
01089     // FlowElem    DOF info
01090     for (size_t i=0; i<D_flow.NodeEssential.size(); i++)
01091     {
01092         D.NodeEssential.push_back(D_flow.NodeEssential[i]);
01093         D.NodeNatural  .push_back(D_flow.NodeNatural  [i]);
01094     }
01095     for (size_t i=0; i<D_flow.FaceEssential.size(); i++)
01096     {
01097         D.FaceEssential.push_back(D_flow.FaceEssential[i]);
01098         D.FaceNatural  .push_back(D_flow.FaceNatural  [i]);
01099     }
01100 
01101     // Insert into DOFInfoMap
01102     DOFInfoMap["Coupled"] = D;
01103 
01104     return 0;
01105 } // }}}
01106 
01107 // Execute the autoregistration
01108 int __CoupledElemDOFInfo_dummy_int  = CoupledDOFInfoRegister();
01109 
01110 }; // namespace FEM
01111 
01112 #endif // FEM_COUPLED_H
01113 
01114 // vim:fdm=marker
Generated on Wed Jan 24 15:56:25 2007 for MechSys by  doxygen 1.4.7