MechSys  1.0
Computing library for simulations in continuum and discrete mechanics
/home/dorival/mechsys/lib/fem/element.h
Go to the documentation of this file.
00001 /************************************************************************
00002  * MechSys - Open Library for Mechanical Systems                        *
00003  * Copyright (C) 2005 Dorival M. Pedroso, Raul Durand                   *
00004  * Copyright (C) 2009 Sergio Galindo                                    *
00005  *                                                                      *
00006  * This program is free software: you can redistribute it and/or modify *
00007  * it under the terms of the GNU General Public License as published by *
00008  * the Free Software Foundation, either version 3 of the License, or    *
00009  * any later version.                                                   *
00010  *                                                                      *
00011  * This program is distributed in the hope that it will be useful,      *
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of       *
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the         *
00014  * GNU General Public License for more details.                         *
00015  *                                                                      *
00016  * You should have received a copy of the GNU General Public License    *
00017  * along with this program. If not, see <http://www.gnu.org/licenses/>  *
00018  ************************************************************************/
00019 
00020 #ifndef MECHSYS_FEM_ELEMENT_H
00021 #define MECHSYS_FEM_ELEMENT_H
00022 
00023 // MechSys
00024 #include <mechsys/geomtype.h>
00025 #include <mechsys/mesh/mesh.h>
00026 #include <mechsys/fem/node.h>
00027 #include <mechsys/fem/geomelem.h>
00028 #include <mechsys/models/model.h>
00029 #include <mechsys/util/maps.h>
00030 #include <mechsys/util/fatal.h>
00031 #include <mechsys/linalg/matvec.h>
00032 
00033 namespace FEM
00034 {
00035 
00036 class Element; 
00037 
00038 class MPyPrms 
00039 {
00040 public:
00041     mutable double          SF;          
00042     mutable double          MaxDist;     
00043     double                  PctMaxDist;  
00044     bool                    AutoLimits;  
00045     bool                    PNG;         
00046     char const            * Extra;       
00047     bool                    OnlyBeams;   
00048     mutable Element const * EleMmin;     
00049     mutable Element const * EleMmax;     
00050     mutable Element const * EleNmin;     
00051     mutable Element const * EleNmax;     
00052     mutable Element const * EleVmin;     
00053     mutable Element const * EleVmax;     
00054     mutable double          Mmin;        
00055     mutable double          Mmax;        
00056     mutable double          Nmin;        
00057     mutable double          Nmax;        
00058     mutable double          Vmin;        
00059     mutable double          Vmax;        
00060     mutable double          rMmin;       
00061     mutable double          rMmax;       
00062     bool                    DrawIPs;     
00063     bool                    FindMLimits; 
00064     size_t                  NDiv;        
00065     bool                    WithTxt;     
00066     bool                    OnlyTxtLim;  
00067     bool                    DrawN;       
00068     bool                    DrawV;       
00069     size_t                  TxtSz;       
00070     bool                    WithReac;    
00071     Array<int>              ReacNodes;   
00072     MPyPrms ()                           
00073     {
00074         SF          = 1.0;
00075         MaxDist     = 1.0;
00076         PctMaxDist  = 0.1;
00077         AutoLimits  = true;
00078         PNG         = false;
00079         Extra       = NULL;
00080         OnlyBeams   = false;
00081         EleMmin     = NULL;
00082         EleMmax     = NULL;
00083         EleNmin     = NULL;
00084         EleNmax     = NULL;
00085         EleVmin     = NULL;
00086         EleVmax     = NULL;
00087         Mmin        = 0.0;
00088         Mmax        = 0.0;
00089         Nmin        = 0.0;
00090         Nmax        = 0.0;
00091         Vmin        = 0.0;
00092         Vmax        = 0.0;
00093         rMmin       = 0.0;
00094         rMmax       = 0.0;
00095         DrawIPs     = true;
00096         FindMLimits = true;
00097         NDiv        = 10;
00098         WithTxt     = true;
00099         OnlyTxtLim  = false;
00100         DrawN       = false;
00101         DrawV       = false;
00102         TxtSz       = 8;
00103         WithReac    = true;
00104     }
00105 };
00106 
00107 
00108 class Element
00109 {
00110 public:
00111     // Constructor
00112     Element (int                  NDim,   
00113              Mesh::Cell   const & Cell,   
00114              Model        const * Mdl,    
00115              Model        const * XMdl,   
00116              SDPair       const & Prp,    
00117              SDPair       const & Ini,    
00118              Array<Node*> const & Nodes); 
00119 
00120     // Destructor
00121     virtual ~Element ();
00122 
00123     // Methods
00124     virtual void IncNLocDOF   (size_t & NEq)                                      const {} 
00125     virtual void BackupState  ()                                                  const;   
00126     virtual void RestoreState ()                                                  const;   
00127     virtual void SetBCs       (size_t IdxEdgeOrFace, SDPair const & BCs, BCFuncs * BCF) {} 
00128     virtual void AddToF       (double Time, Vec_t & F)                            const {} 
00129     virtual void ClrBCs       ()                                                        {} 
00130     virtual void GetLoc       ()                                                  const { throw new Fatal("Element::GetLoc: (no arguments) Method not implemented for this element"); } 
00131     virtual void GetLoc       (Array<size_t> & Loc)                               const { throw new Fatal("Element::GetLoc: Method not implemented for this element"); } 
00132     virtual void ElemEqs      (double Time, Vec_t const & V, Vec_t const & Ww, Vec_t const & Pw) const { throw new Fatal("Element::GetLoc: Method not implemented for this element"); } 
00133     virtual void CalcK        (Mat_t & K)                                         const { throw new Fatal("Element::CalcK: Method not implemented for this element"); }
00134     virtual void CalcM        (Mat_t & M)                                         const { throw new Fatal("Element::CalcM: Method not implemented for this element"); }
00135     virtual void CalcC        (Mat_t & C)                                         const { throw new Fatal("Element::CalcC: Method not implement for this element"); }
00136     virtual void CalcK        (Vec_t const & U, double Alpha, double dt, Mat_t & KK, Vec_t & dF) const { throw new Fatal("Element::CalcK (with U, Alpha, Dt => HydroMech): Method not implemented for this element"); }
00137     virtual void CalcKCM      (Mat_t & KK, Mat_t & CC, Mat_t & MM)                const { throw new Fatal("Element::CalcKCM: Method not implemented for this element"); }
00138     virtual void UpdateState  (Vec_t const & dU, Vec_t * F_int=NULL)              const {}
00139     virtual void SetFint      (Vec_t * Fint=NULL)                                 const {} 
00140     virtual void StateKeys    (Array<String> & Keys)                              const {} 
00141     virtual void StateAtIP    (SDPair & KeysVals, int IdxIP)                      const {} 
00142     virtual void StateAtIPs   (Array<SDPair> & Results)                           const;   
00143     virtual void StateAtNodes (Array<SDPair> & Results)                           const;   
00144     virtual void Draw         (std::ostream & os, MPyPrms const & Prms)           const;   
00145 
00146     // Methods for the Runge-Kutta method
00147     virtual size_t NIVs       ()                                                            const { return 0; } 
00148     virtual double GetIV      (size_t i)                                                    const { return 0; } 
00149     virtual void   SetIV      (size_t i, double Val)                                              {}            
00150     virtual void   CalcIVRate (double Time, Vec_t const & U, Vec_t const & V, Vec_t & Rate) const {}            
00151     virtual void   CorrectIVs ()                                                                  {}            
00152 
00153     virtual double Update  (size_t Idx, Vec_t const & V_g, Vec_t const & dPwdt_g, double dt) { throw new Fatal("Element::Update: Method not implemented in this element"); }  
00154     virtual void   Restore ()                                                                { throw new Fatal("Element::Restore: Method not implemented in this element"); } 
00155 
00156     // Methods that depend on GE
00157     void CoordMatrix   (Mat_t & C)                 const; 
00158     void FCoordMatrix  (size_t IdxFace, Mat_t & C) const; 
00159     void ShapeMatrix   (Mat_t & M)                 const; 
00160     void CalcShape     (Mat_t const & C,  IntegPoint const & IP,  double & detJ, double & Coef) const; 
00161     void CalcFaceShape (Mat_t const & FC, IntegPoint const & FIP, double & detJ, double & Coef, double h=1.0) const; 
00162     void CalcFaceShape (Mat_t const & FC, IntegPoint const & FIP, Mat_t & J, double & detJ, double & Coef, double h=1.0) const; 
00163     void CoordsOfIP    (size_t IdxIP, Vec_t & X)   const; 
00164 
00165     // Data
00166     int                NDim;   
00167     Mesh::Cell const & Cell;   
00168     Model      const * Mdl;    
00169     Model      const * XMdl;   
00170     GeomElem         * GE;     
00171     bool               Active; 
00172     GeomType           GTy;    
00173     Array<Node*>       Con;    
00174     Array<State*>      Sta;    
00175     bool               IsUWP;  
00176 };
00177 
00178 
00180 
00181 
00182 inline Element::Element (int TheNDim, Mesh::Cell const & TheCell, Model const * TheMdl, Model const * TheXMdl, SDPair const & Prp, SDPair const & Ini, Array<Node*> const & Nodes)
00183     : NDim(TheNDim), Cell(TheCell), Mdl(TheMdl), XMdl(TheXMdl), GE(NULL), Active(Prp.HasKey("active") ? Prp("active") : true),
00184       GTy(SDPairToGType(Prp,(NDim==3?"d3d":"d2d"))), IsUWP(false)
00185 {
00186     // connectivity
00187     Con.Resize (Nodes.Size());
00188     for (size_t i=0; i<Nodes.Size(); ++i)
00189     {
00190         Con[i] = Nodes[i];
00191         if (Active) Con[i]->NShares++;
00192     }
00193 
00194     // geometry element
00195     if (Prp.HasKey("geom"))
00196     {
00197         String geom_name;
00198         GEOM.Val2Key (Prp("geom"), geom_name);
00199         GE = AllocGeomElem (geom_name, NDim);
00200 
00201         // set number of integration points
00202         if (Prp.HasKey("nip")) GE->SetIPs (static_cast<int>(Prp("nip")));
00203 
00204         // check number of nodes
00205         if (Con.Size()!=GE->NN) throw new Fatal("Element::Element: Number of vertices (%d) of an element (%s) in mesh is incorrect. It must be equal to %d for %s", Con.Size(),geom_name.CStr(),GE->NN,geom_name.CStr());
00206     }
00207 
00208     // check model
00209     if (Mdl!=NULL) { if (GTy!=Mdl->GTy) throw new Fatal("Element::Element: Element geometry type (%s) must be equal to Model geometry type (%s)", GTypeToStr(GTy).CStr(), GTypeToStr(Mdl->GTy).CStr()); }
00210 }
00211 
00212 inline Element::~Element ()
00213 {
00214     if (GE!=NULL) delete GE;
00215     for (size_t i=0; i<Sta.Size(); ++i) { if (Sta[i]!=NULL) delete Sta[i]; }
00216 }
00217 
00218 inline void Element::BackupState () const
00219 {
00220     for (size_t i=0; i<Sta.Size(); ++i) Sta[i]->Backup();
00221 }
00222 
00223 inline void Element::RestoreState () const 
00224 {
00225     for (size_t i=0; i<Sta.Size(); ++i) Sta[i]->Restore();
00226 }
00227 
00228 inline void Element::CoordMatrix (Mat_t & C) const
00229 {
00230     if (GE==NULL) throw new Fatal("Element::CoordMatrix: This method works only when GE (geometry element) is not NULL");
00231 
00232     C.change_dim (GE->NN, NDim);
00233     for (size_t i=0; i<GE->NN; ++i)
00234     for (int    j=0; j<NDim;   ++j)
00235         C(i,j) = Con[i]->Vert.C[j];
00236 }
00237 
00238 inline void Element::FCoordMatrix (size_t IdxFace, Mat_t & C) const
00239 {
00240     if (GE==NULL) throw new Fatal("Element::FCoordMatrix: This method works only when GE (geometry element) is not NULL");
00241 
00242     C.change_dim (GE->NFN, NDim);
00243     for (size_t i=0; i<GE->NFN; ++i)
00244     for (int    j=0; j<NDim;    ++j)
00245         C(i,j) = Con[GE->FNode(IdxFace,i)]->Vert.C[j];
00246 }
00247 
00248 inline void Element::ShapeMatrix (Mat_t & M) const
00249 {
00250     if (GE==NULL) throw new Fatal("Element::ShapeMatrix: This method works only when GE (geometry element) is not NULL");
00251 
00252     /* Shape functions matrix:
00253               _                                              _
00254              |   N11      N12      N13      ...  N1(NN)       |
00255              |   N21      N22      N23      ...  N2(NN)       |
00256          M = |   N31      N32      N33      ...  N3(NN)       |
00257              |          ...............................       |
00258              |_  N(NIP)1  N(NIP)2  N(NIP)3  ...  N(NIP)(NN)  _| (NIP x NN)
00259     */
00260     M.change_dim (GE->NIP, GE->NN);
00261     for (size_t i=0; i<GE->NIP; ++i)
00262     {
00263         GE->Shape (GE->IPs[i].r, GE->IPs[i].s, GE->IPs[i].t);
00264         for (size_t j=0; j<GE->NN; ++j) M(i,j) = GE->N(j);
00265     }
00266 }
00267 
00268 inline void Element::CalcShape (Mat_t const & C, IntegPoint const & IP, double & detJ, double & Coef) const
00269 {
00270     if (GE==NULL) throw new Fatal("Element::ShapeMatrix: This method works only when GE (geometry element) is not NULL");
00271 
00272     // geometric data
00273     GE->Shape  (IP.r, IP.s, IP.t);
00274     GE->Derivs (IP.r, IP.s, IP.t);
00275 
00276     // Jacobian and its determinant
00277     Mat_t J(GE->dNdR * C); // J = dNdR * C
00278     detJ = Det(J);
00279 
00280     // coefficient used during integration
00281     Coef = detJ*IP.w;
00282 }
00283 
00284 inline void Element::CalcFaceShape (Mat_t const & FC, IntegPoint const & FIP, double & detJ, double & Coef, double h) const
00285 {
00286     if (GE==NULL) throw new Fatal("Element::ShapeMatrix: This method works only when GE (geometry element) is not NULL");
00287 
00288     // geometric data
00289     GE->FaceShape  (FIP.r, FIP.s);
00290     GE->FaceDerivs (FIP.r, FIP.s);
00291 
00292     // face/edge Jacobian and its determinant
00293     Mat_t J(GE->FdNdR * FC);
00294     detJ = Det(J);
00295 
00296     // coefficient used during integration
00297     Coef = h*detJ*FIP.w;
00298 
00299     // correct Coef for axisymmetric problems
00300     if (GTy==axs_t)
00301     {
00302         // calculate radius=x at this FIP
00303         double radius = 0.0;
00304         for (size_t i=0; i<GE->NFN; ++i) radius += GE->FN(i)*FC(i,0);
00305         Coef *= radius;
00306     }
00307 }
00308 
00309 inline void Element::CalcFaceShape (Mat_t const & FC, IntegPoint const & FIP, Mat_t & J, double & detJ, double & Coef, double h) const
00310 {
00311     if (GE==NULL) throw new Fatal("Element::ShapeMatrix: This method works only when GE (geometry element) is not NULL");
00312 
00313     // geometric data
00314     GE->FaceShape  (FIP.r, FIP.s);
00315     GE->FaceDerivs (FIP.r, FIP.s);
00316 
00317     // face/edge Jacobian and its determinant
00318     //J.change_dim (NDim-1, NDim);
00319     J    = GE->FdNdR * FC;
00320     detJ = Det(J);
00321 
00322     // coefficient used during integration
00323     Coef = h*detJ*FIP.w;
00324 
00325     // correct Coef for axisymmetric problems
00326     if (GTy==axs_t)
00327     {
00328         // calculate radius=x at this FIP
00329         double radius = 0.0;
00330         for (size_t i=0; i<GE->NFN; ++i) radius += GE->FN(i)*FC(i,0);
00331         Coef *= radius;
00332     }
00333 }
00334 
00335 inline void Element::StateAtIPs (Array<SDPair> & Results) const
00336 {
00337     if (GE==NULL) throw new Fatal("Element::StateAtIPs: This method works only when GE (geometry element) is not NULL");
00338 
00339     // one set of results per IP
00340     Results.Resize (GE->NIP);
00341     for (size_t i=0; i<GE->NIP; ++i) StateAtIP (Results[i], i);
00342 }
00343 
00344 inline void Element::StateAtNodes (Array<SDPair> & Results) const
00345 {
00346     if (GE==NULL) throw new Fatal("Element::StateAtNodes: This method works only when GE (geometry element) is not NULL");
00347 
00348     // resize results array (one set of results per node)
00349     Results.Resize (GE->NN);
00350 
00351     // matrix of extrapolation
00352     Mat_t E;
00353     if (GE->NIP < GE->NN)
00354     {
00355         // matrix with natural coordinates of nodes
00356         Mat_t Cnat;
00357         GE->NatCoords (Cnat); // size = GE->NN x NDim+1
00358 
00359         // auxiliar matrices
00360         Mat_t M(GE->NIP, GE->NN); // shape func matrix
00361         Mat_t C(GE->NIP, NDim+1); // matrix with coordinates of IPs
00362         for (size_t i=0; i<GE->NIP; ++i)
00363         {
00364             GE->Shape (GE->IPs[i].r, GE->IPs[i].s, GE->IPs[i].t);
00365             for (size_t j=0; j<GE->NN; ++j) M(i,j) = GE->N(j);
00366             C(i,0) = GE->IPs[i].r;
00367             C(i,1) = GE->IPs[i].s;
00368             if (NDim==2) C(i,2) = 1.0;
00369             else
00370             {
00371                 C(i,2) = GE->IPs[i].t;
00372                 C(i,3) = 1.0;
00373             }
00374         }
00375 
00376         // matrix of extrapolation
00377         Mat_t I, Mi, Ci;
00378         Identity (GE->NIP, I);
00379         Inv (M, Mi);
00380         Inv (C, Ci);
00381         Mat_t tmp(I - C*Ci);
00382         E = Mi*tmp + Cnat*Ci;
00383     }
00384     else
00385     {
00386         // shape func matrix
00387         Mat_t M;
00388         ShapeMatrix (M);
00389         Inv (M, E);
00390     }
00391 
00392     // keys
00393     Array<String> keys;
00394     StateKeys (keys);
00395 
00396     // state at IPs
00397     Array<SDPair> state_at_IPs;
00398     StateAtIPs (state_at_IPs);
00399 
00400     // extrapolate
00401     Vec_t val_at_IPs (GE->NIP); // values at IPs
00402     Vec_t val_at_Nods(GE->NN);  // values at nodes
00403     for (size_t k=0; k<keys.Size(); ++k)
00404     {
00405         // gather values from IPs
00406         for (size_t i=0; i<GE->NIP; ++i) val_at_IPs(i) = state_at_IPs[i](keys[k]);
00407 
00408         // extrapolate to nodes
00409         val_at_Nods = E * val_at_IPs;
00410 
00411         // scatter values to nodes
00412         for (size_t i=0; i<GE->NN; ++i) Results[i].Set (keys[k].CStr(), val_at_Nods(i));
00413     }
00414 }
00415 
00416 inline void Element::CoordsOfIP (size_t IdxIP, Vec_t & X) const
00417 {
00418     if (GE==NULL) throw new Fatal("Element::CoordsOfIP: This method works only when GE (geometry element) is not NULL");
00419 
00420     //std::cout << "\n" << GE->IPs[IdxIP].r << "  " << GE->IPs[IdxIP].s << "\n\n";
00421     GE->Shape (GE->IPs[IdxIP].r, GE->IPs[IdxIP].s, GE->IPs[IdxIP].t);
00422     X.change_dim (NDim);
00423     set_to_zero  (X);
00424     for (size_t i=0; i<GE->NN; ++i)
00425     for (int    j=0; j<NDim;   ++j)
00426         X(j) += GE->N(i)*Con[i]->Vert.C[j];
00427 }
00428 
00429 inline void Element::Draw (std::ostream & os, MPyPrms const & Prms) const
00430 {
00431     if (GE==NULL) throw new Fatal("Element::CoordsIP: This method works only when GE (geometry element) is not NULL");
00432 
00433     // draw shape
00434     os << "dat.append((PH.MOVETO, (" << Con[0]->Vert.C[0] << "," << Con[0]->Vert.C[1] << ")))\n";
00435     if (GE->NN<=4)
00436     {
00437         for (size_t i=1; i<GE->NN; ++i) os << "dat.append((PH.LINETO,    (" << Con[i]->Vert.C[0] << "," << Con[i]->Vert.C[1] << ")))\n";
00438         os << "dat.append((PH.CLOSEPOLY, (" << Con[0]->Vert.C[0] << "," << Con[0]->Vert.C[1] << ")))\n";
00439     }
00440     else if (GE->NN==6)
00441     {
00442         os << "dat.append((PH.LINETO,    (" << Con[3]->Vert.C(0) << "," << Con[3]->Vert.C(1) << ")))\n";
00443         os << "dat.append((PH.LINETO,    (" << Con[1]->Vert.C(0) << "," << Con[1]->Vert.C(1) << ")))\n";
00444         os << "dat.append((PH.LINETO,    (" << Con[4]->Vert.C(0) << "," << Con[4]->Vert.C(1) << ")))\n";
00445         os << "dat.append((PH.LINETO,    (" << Con[2]->Vert.C(0) << "," << Con[2]->Vert.C(1) << ")))\n";
00446         os << "dat.append((PH.LINETO,    (" << Con[5]->Vert.C(0) << "," << Con[5]->Vert.C(1) << ")))\n";
00447         os << "dat.append((PH.CLOSEPOLY, (" << Con[0]->Vert.C(0) << "," << Con[0]->Vert.C(1) << ")))\n";
00448     }
00449     else if (GE->NN==8)
00450     {
00451         os << "dat.append((PH.LINETO,    (" << Con[4]->Vert.C(0) << "," << Con[4]->Vert.C(1) << ")))\n";
00452         os << "dat.append((PH.LINETO,    (" << Con[1]->Vert.C(0) << "," << Con[1]->Vert.C(1) << ")))\n";
00453         os << "dat.append((PH.LINETO,    (" << Con[5]->Vert.C(0) << "," << Con[5]->Vert.C(1) << ")))\n";
00454         os << "dat.append((PH.LINETO,    (" << Con[2]->Vert.C(0) << "," << Con[2]->Vert.C(1) << ")))\n";
00455         os << "dat.append((PH.LINETO,    (" << Con[6]->Vert.C(0) << "," << Con[6]->Vert.C(1) << ")))\n";
00456         os << "dat.append((PH.LINETO,    (" << Con[3]->Vert.C(0) << "," << Con[3]->Vert.C(1) << ")))\n";
00457         os << "dat.append((PH.LINETO,    (" << Con[7]->Vert.C(0) << "," << Con[7]->Vert.C(1) << ")))\n";
00458         os << "dat.append((PH.CLOSEPOLY, (" << Con[0]->Vert.C(0) << "," << Con[0]->Vert.C(1) << ")))\n";
00459     }
00460 
00461     // model has deviatoric plastic strain ?
00462     bool has_edp = false;
00463     if (Mdl->IvNames.Find("edp")>=0) has_edp = true;
00464 
00465     // draw IPs
00466     if (Prms.DrawIPs)
00467     {
00468         os << "x_ips, y_ips = [], []\n";
00469         os << "x_edp, y_edp = [], []\n";
00470         for (size_t i=0; i<GE->NIP; ++i)
00471         {
00472             Vec_t X;
00473             CoordsOfIP (i,X);
00474             os << "x_ips.append(" << X(0) << ")\n";
00475             os << "y_ips.append(" << X(1) << ")\n";
00476             if (has_edp)
00477             {
00478                 SDPair res;
00479                 StateAtIP (res, i);
00480                 if (res("edp")>0.0)
00481                 {
00482                     os << "x_edp.append(" << X(0) << ")\n";
00483                     os << "y_edp.append(" << X(1) << ")\n";
00484                 }
00485             }
00486         }
00487         os << "plot(x_ips,y_ips,'r*',zorder=1)\n";
00488     }
00489     if (has_edp) os << "plot(x_edp,y_edp,'ko',ms=9,zorder=2)\n";
00490 }
00491 
00492 std::ostream & operator<< (std::ostream & os, Element const & E)
00493 {
00494     os << Util::_4 << E.Cell.ID << " " << Util::_4 << E.Cell.Tag << " ";
00495     os << (E.Active?TERM_GREEN:TERM_RED) << (E.Active?" active":" inactive") << TERM_RST << " ";
00496     os << GTypeToStr(E.GTy) << " ";
00497     if (E.GE!=NULL)  os << E.GE->Name  << " " << "NIP=" << E.GE->NIP << " ";
00498     else             os << "GE=NULL"   << " ";
00499     if (E.Mdl!=NULL) os << E.Mdl->Name << " ";
00500     else             os << "Mdl=NULL";
00501     os << "(";
00502     for (size_t i=0; i<E.Con.Size(); ++i)
00503     {
00504         os << E.Con[i]->Vert.ID;
00505         if (i==E.Con.Size()-1) os << ") ";
00506         else                   os << ",";
00507     }
00508     return os;
00509 }
00510 
00511 
00513 
00514 
00515 typedef Element * (*ElementMakerPtr)(int NDim, Mesh::Cell const & Cell, Model const * Mdl, Model const * XMdl, SDPair const & Prp, SDPair const & Ini, Array<Node*> const & Nodes);
00516 
00517 typedef std::map<String, ElementMakerPtr> ElementFactory_t;
00518 
00519 ElementFactory_t ElementFactory;
00520 
00521 Element * AllocElement(String const & Name, int NDim, Mesh::Cell const & Cell, Model const * Mdl, Model const * XMdl, SDPair const & Prp, SDPair const & Ini, Array<Node*> const & Nodes)
00522 {
00523     ElementFactory_t::iterator it = ElementFactory.find(Name);
00524     if (it==ElementFactory.end()) throw new Fatal("AllocElement: '%s' is not available", Name.CStr());
00525 
00526     Element * ptr = (*it->second)(NDim,Cell,Mdl,XMdl,Prp,Ini,Nodes);
00527 
00528     return ptr;
00529 }
00530 
00531 
00533 
00534 
00535 SDPair PROB; 
00536 
00537 typedef std::map<String, std::pair<String,String> > ElementVarKeys_t;   
00538 typedef std::map<String, Array<String> >            ElementExtraKeys_t; 
00539 
00540 ElementVarKeys_t   ElementVarKeys;   
00541 ElementExtraKeys_t ElementExtraKeys; 
00542 
00543 }; // namespace FEM
00544 
00545 #ifdef USE_BOOST_PYTHON
00546 double PyPROB (BPy::str const & Key) { return FEM::PROB(BPy::extract<char const *>(Key)()); }
00547 #endif
00548 
00549 #endif // MECHSYS_FEM_ELEMENT
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines