00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #ifndef MECHSYS_FEM_STAGESMANAGER_H
00023 #define MECHSYS_FEM_STAGESMANAGER_H
00024
00025 #ifdef HAVE_CONFIG_H
00026 #include "config.h"
00027 #else
00028 #ifndef REAL
00029 #define REAL double
00030 #endif
00031 #endif
00032
00033 #include <algorithm>
00034
00035 #include "util/string.h"
00036 #include "fem/inputdata.h"
00037 #include "fem/filesdata.h"
00038 #include "fem/data.h"
00039
00040 namespace FEM
00041 {
00042
00044
00052 class StagesManager
00053 {
00054 public:
00055
00056 friend std::ostream & operator<< (std::ostream & os, FEM::StagesManager const & SM);
00057
00058 void Initialize (InputData * ID, FilesData * FDat, Data * data);
00059 int nStages () const { return _idat->nSTAGES; }
00060 int CurrentStageNum () const { return _current_stage_num; }
00061 void ActivateStage (int iStage);
00062 void DumpStage (int nStage) {}
00063 void RestoreStage (int nStage) {}
00064 private:
00065
00066 int _current_stage_num;
00067 InputData * _idat;
00068 FilesData * _fdat;
00069 Data * _data;
00070
00071 void _alloc_and_fill_nodal_coords();
00072 void _alloc_and_init_elements();
00073 void _alloc_and_init_embededs();
00074 void _clear_bry_conditions();
00075 void _update_nodes_given_face_bry(int iStage);
00076 void _update_nodes_given_node_bry(int iStage);
00077 void _set_element_attributes(int iStage);
00078 };
00079
00080
00082
00083
00084 inline void StagesManager::Initialize(InputData * ID, FilesData * FDat, Data * data)
00085 {
00086 _current_stage_num = -1;
00087 _idat = ID;
00088 _fdat = FDat;
00089 _data = data;
00090 _alloc_and_fill_nodal_coords();
00091 _alloc_and_init_elements();
00092 }
00093
00094 inline void StagesManager::_alloc_and_fill_nodal_coords()
00095 {
00096
00097 _data->Nodes.resize(_fdat->Nodes.size());
00098
00099
00100 for (size_t i=0; i<_data->Nodes.size(); ++i)
00101 {
00102
00103 _data->Nodes[i].Initialize( i,
00104 _fdat->Nodes[i].X,
00105 _fdat->Nodes[i].Y,
00106 _fdat->Nodes[i].Z);
00107 }
00108 }
00109
00110 inline void StagesManager::_alloc_and_init_elements()
00111 {
00112
00113 _data->Elements.resize(_fdat->Elements.size());
00114
00115
00116 FilesData::Stage const & S = _fdat->Stages[0];
00117
00118
00119 for (size_t i=0; i<_data->Elements.size(); ++i)
00120 {
00121
00122 FilesData::Element const & E = _fdat->Elements[i];
00123
00124
00125
00126 _data->Elements[i] = AllocElement(S.Attributes[E.IdxAtt].EleName);
00127 _data->Elements[i]->Number = i+1;
00128
00129
00130 for (size_t j=0; j<E.Connects.size(); ++j)
00131 {
00132
00133
00134
00135 _data->Elements[i]->SetNode(static_cast<int>(j), &_data->Nodes[E.Connects[j]-1]);
00136 _data->Nodes[E.Connects[j]-1].Refs()++;
00137 }
00138 }
00139
00140
00141 if (_fdat->EmbElements.size()>0) _alloc_and_init_embededs();
00142
00143 }
00144
00145 inline void StagesManager::_alloc_and_init_embededs()
00146 {
00147
00148
00149
00150 FilesData::Stage const & S = _fdat->Stages[0];
00151 int n_elems = _fdat->Elements.size();
00152 int n_embs = _fdat->EmbElements.size();
00153 REAL tiny = 1E-4;
00154
00155 for (int i_emb=0; i_emb<n_embs; i_emb++)
00156 {
00157 FilesData::Element & great_embed = _fdat->EmbElements[i_emb];
00158 great_embed.IniVals.resize(1);
00159 great_embed.IniVals[0].resize(1);
00160 great_embed.IniVals[0][0]=0.0;
00161 FilesData::Node & init_node = _fdat->Nodes[great_embed.Connects[0]-1];
00162 FilesData::Node & final_node = _fdat->Nodes[great_embed.Connects[1]-1];
00163 REAL x1 = init_node.X;
00164 REAL y1 = init_node.Y;
00165 REAL z1 = init_node.Z;
00166 REAL x2 = final_node.X;
00167 REAL y2 = final_node.Y;
00168 REAL z2 = final_node.Z;
00169
00170 REAL length = sqrt(pow(x2-x1,2)+pow(y2-y1,2)+pow(z2-z1,2));
00171 REAL tinylength=0.01*length;
00172 REAL l = (x2-x1)/length;
00173 REAL m = (y2-y1)/length;
00174 REAL n = (z2-z1)/length;
00175 int init_elem = 0;
00176 int prev_elem = 0;
00177 int actual_elem = 0;
00178 int next_elem = 0;
00179 int final_elem = 0;
00180
00181 for (int j=0; j<n_elems; j++)
00182 if (_data->Elements[j]->IsInside(x1+tinylength*l, y1+tinylength*m, z1+tinylength*n))
00183 {
00184 init_elem = j; break;
00185 }
00186
00187 for (int j=0; j<n_elems; j++)
00188 if (_data->Elements[j]->IsInside(x2-tinylength*l, y2-tinylength*m, z2-tinylength*n))
00189 {
00190 final_elem = j; break;
00191 }
00192 REAL xp=x1;
00193 REAL yp=y1;
00194 REAL zp=z1;
00195 REAL x=xp;
00196 REAL y=yp;
00197 REAL z=zp;
00198 actual_elem = prev_elem = init_elem;
00199 next_elem = 0;
00200 REAL step = length;
00201 bool final = false;
00202
00203 if (init_elem == final_elem)
00204 {
00205 x=x2; y=y2; z=z2; final=true;
00206 }
00207
00208
00209 do{
00210 if (!final) {
00211 step *= 0.501;
00212 x += step*l;
00213 y += step*m;
00214 z += step*n;
00215 for (int j=0; j<n_elems; j++)
00216 if (_data->Elements[j]->IsInside(x, y, z))
00217 {
00218 actual_elem = j; break;
00219 }
00220 }
00221 if (prev_elem == actual_elem)
00222 {
00223 REAL bf;
00224 if (!final)
00225 {
00226 REAL r, s, t;
00227 _data->Elements[next_elem]->InverseMap(x,y,z,r,s,t);
00228 bf = fabs(_data->Elements[next_elem]->BoundDistance(r,s,t));
00229 }
00230 if ( final || bf <= tiny )
00231 {
00232
00233 Element * elem = AllocElement(S.Attributes[great_embed.IdxAtt].EleName);
00234 Embedded * embedded = dynamic_cast<Embedded*>(elem);
00235
00236
00237 Array<Node*> nodes;
00238 size_t n_emb_nodes = elem->nNodes();
00239 nodes.resize(n_emb_nodes);
00240
00241
00242 size_t n_div = n_emb_nodes - 1;
00243 for (size_t j=0; j<n_emb_nodes; ++j)
00244 {
00245 nodes[j] = new Node;
00246 REAL _x = xp + j*(x-xp)/n_div;
00247 REAL _y = yp + j*(y-yp)/n_div;
00248 REAL _z = zp + j*(z-zp)/n_div;
00249 nodes[j]->Initialize(0, _x, _y, _z);
00250 }
00251
00252
00253 embedded->Initialize(nodes, _data->Elements[actual_elem]);
00254
00255 _data->Elements.push_back(embedded);
00256 _fdat->Elements.push_back(_fdat->EmbElements[i_emb]);
00257
00258 if (final) break;
00259 else{
00260 xp=x; yp=y; zp=z;
00261 if (next_elem==final_elem)
00262 {
00263 x=x2; y=y2; z=z2;
00264 actual_elem = prev_elem = next_elem;
00265 final = true;
00266 }
00267 else
00268 {
00269 prev_elem = next_elem;
00270 step = sqrt(pow(x-x2,2)+pow(y-y2,2)+pow(z-z2,2));
00271 }
00272 }
00273 }else{
00274 step=fabs(step);
00275 }
00276 }else{
00277 step=-fabs(step);
00278 next_elem = actual_elem;
00279 }
00280 }while (true);
00281 }
00282 }
00283
00284 inline void StagesManager::_clear_bry_conditions()
00285 {
00286 for (size_t i=0; i<_data->Nodes.size(); ++i)
00287 _data->Nodes[i].ClearBryValues();
00288 }
00289
00290 inline void StagesManager::_update_nodes_given_face_bry(int iStage)
00291 {
00292
00293 FilesData::Stage const & S = _fdat->Stages[iStage];
00294
00295
00296 for (size_t i=0; i<_fdat->Faces.size(); ++i)
00297 {
00298
00299 int idx_bry = _fdat->Faces[i].IdxBry;
00300 if (idx_bry!=-1)
00301 {
00302
00303 FilesData::BryCond const & BC = S.FaceBrys[idx_bry];
00304
00305
00306 int n_face_nodes = _fdat->Faces[i].Connects.size();
00307
00308
00309 int k_element = _fdat->Faces[i].OwnerElement-1;
00310
00311
00312 FEM::Element * const ptrE = _data->Elements[k_element];
00313
00314
00315 Array<FEM::Node*> a_ptr_face_nodes;
00316 a_ptr_face_nodes.resize(n_face_nodes);
00317 for (int j=0; j<n_face_nodes; ++j)
00318 {
00319
00320 int j_node = _fdat->Faces[i].Connects[j]-1;
00321 a_ptr_face_nodes[j] = &_data->Nodes[j_node];
00322 }
00323
00324
00325 for (size_t j=0; j<BC.Names.size(); ++j)
00326 {
00327
00328 if (ptrE->IsEssential(BC.Names[j]))
00329 {
00330
00331 for (int m=0; m<n_face_nodes; ++m)
00332 {
00333
00334 FEM::Node::DOFVarsStruct & dof_var = a_ptr_face_nodes[m]->DOFVar(BC.Names[j]);
00335
00336
00337 dof_var.EssentialBry = BC.Values[j];
00338 dof_var.IsEssenPresc = true;
00339 }
00340 }
00341 else
00342 {
00343
00344 String str_dof_name;
00345
00346
00347 LinAlg::Vector<REAL> a_nodal_values;
00348
00349
00350 ptrE->CalcFaceNodalValues(BC.Names[j] ,
00351 BC.Values[j] ,
00352 a_ptr_face_nodes,
00353 str_dof_name ,
00354 a_nodal_values );
00355
00356
00357 for (int m=0; m<n_face_nodes; ++m)
00358 {
00359
00360 FEM::Node::DOFVarsStruct & dof_var = a_ptr_face_nodes[m]->DOFVar(str_dof_name);
00361
00362
00363 dof_var.NaturalBry += a_nodal_values(m);
00364 dof_var.IsEssenPresc = false;
00365 }
00366 }
00367 }
00368 }
00369
00370
00371 }
00372 }
00373
00374 inline void StagesManager::_update_nodes_given_node_bry(int iStage)
00375 {
00376
00377 FilesData::Stage const & S = _fdat->Stages[iStage];
00378
00379
00380 for (size_t i=0; i<_fdat->Nodes.size(); ++i)
00381 {
00382
00383 int idx_bry = _fdat->Nodes[i].IdxBry;
00384 if (idx_bry!=-1 && _data->Nodes[i].Refs()>0)
00385 {
00386
00387 FilesData::BryCond const & BC = S.NodeBrys[idx_bry];
00388
00389
00390 for (size_t j=0; j<BC.Names.size(); ++j)
00391 {
00392
00393 FEM::Node::DOFVarsStruct & dof_var = _data->Nodes[i].DOFVar(BC.Names[j]);
00394 if (_data->Nodes[i].IsEssential(BC.Names[j]))
00395 {
00396 dof_var.EssentialBry = BC.Values[j];
00397 dof_var.IsEssenPresc = true;
00398 }
00399 else
00400 {
00401 dof_var.NaturalBry += BC.Values[j];
00402 dof_var.IsEssenPresc = false;
00403 }
00404 }
00405 }
00406
00407
00408 }
00409 }
00410
00411 inline void StagesManager::_set_element_attributes(int iStage)
00412 {
00413
00414 FilesData::Stage const & S = _fdat->Stages[iStage];
00415
00416
00417 for (size_t i=0; i<_fdat->Elements.size(); ++i)
00418 {
00419
00420 int idx_att = _fdat->Elements[i].IdxAtt;
00421
00422
00423 FilesData::Attribute const & A = S.Attributes[idx_att];
00424
00425
00426 if (A.IsActive) _data->Elements[i]->Activate();
00427 else _data->Elements[i]->Deactivate();
00428
00429
00430 _data->Elements[i]->SetProperties(A.Properties);
00431
00432
00433 _data->Elements[i]->ReAllocateModel(A.ModelName, A.ModelPrms, _fdat->Elements[i].IniVals);
00434 }
00435 }
00436
00437 inline void StagesManager::ActivateStage(int iStage)
00438 {
00439 _current_stage_num = iStage;
00440 _clear_bry_conditions();
00441 _update_nodes_given_face_bry(iStage);
00442 _update_nodes_given_node_bry(iStage);
00443 _set_element_attributes(iStage);
00444
00445 #ifdef DO_DEBUG
00446 std::cout << "Stage " << iStage << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
00447 std::cout << (*_data) << std::endl;
00448 #endif
00449
00450 }
00451
00452 };
00453
00454 #endif // MECHSYS_FEM_STAGESMANAGER_H
00455
00456