![]() |
MechSys
1.0
Computing library for simulations in continuum and discrete mechanics
|
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