v0.13.1
homogenisation_periodic_atom_test.cpp
Go to the documentation of this file.
1/* This file is part of MoFEM.
2 * MoFEM is free software: you can redistribute it and/or modify it under
3 * the terms of the GNU Lesser General Public License as published by the
4 * Free Software Foundation, either version 3 of the License, or (at your
5 * option) any later version.
6 *
7 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10 * License for more details.
11 *
12 * You should have received a copy of the GNU Lesser General Public
13 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14
15
16#include <MoFEM.hpp>
17using namespace MoFEM;
18
20#include <petsctime.h>
21
22// #include <FEMethod_LowLevelStudent.hpp>
23// #include <FEMethod_UpLevelStudent.hpp>
24
25#include <adolc/adolc.h>
27#include <Hooke.hpp>
28
29#include <boost/shared_ptr.hpp>
30#include <boost/numeric/ublas/vector_proxy.hpp>
31#include <boost/iostreams/tee.hpp>
32#include <boost/iostreams/stream.hpp>
33#include <fstream>
34#include <iostream>
35
37#include <TimeForceScale.hpp>
38#include <VolumeCalculation.hpp>
43
44
45#include <boost/ptr_container/ptr_map.hpp>
46#include <PostProcOnRefMesh.hpp>
47#include <PostProcStresses.hpp>
48
49
50
51
52
53
54static char help[] = "...\n\n";
55
56
57//=================================================================================================================================
58//Define class and multindex container to store data for traiangles on the boundary of the RVE (it cannot be defined within main)
59//=================================================================================================================================
60
63 const EntityHandle Tri_Hand;
64 Face_CenPos_Handle(double _xcoord, double _ycoord, double _zcoord, const EntityHandle _Tri_Hand):xcoord(_xcoord),
65 ycoord(_ycoord), zcoord(_zcoord), Tri_Hand(_Tri_Hand) {}
66};
67
68struct xcoord_tag {}; //tags to used in the multindex container
69struct ycoord_tag {};
70struct zcoord_tag {};
71struct Tri_Hand_tag {};
73
74typedef multi_index_container<
76indexed_by<
77ordered_non_unique<
78tag<xcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord> >,
79
80ordered_non_unique<
81tag<ycoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord> >,
82
83ordered_non_unique<
84tag<zcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> >,
85
86ordered_unique<
87tag<Tri_Hand_tag>, member<Face_CenPos_Handle,const EntityHandle,&Face_CenPos_Handle::Tri_Hand> >,
88
89ordered_unique<
90tag<Composite_xyzcoord>,
91composite_key<
93member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord>,
94member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord>,
95member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> > >
97
98//============================================================================================
99//============================================================================================
100
101
102#define RND_EPS 1e-6
103
104//Rounding
105double roundn(double n)
106{
107 //break n into fractional part (fract) and integral part (intp)
108 double fract, intp;
109 fract = modf(n,&intp);
110
111// //round up
112// if (fract>=.5)
113// {
114// n*=10;
115// ceil(n);
116// n/=10;
117// }
118// //round down
119// if (fract<=.5)
120// {
121// n*=10;
122// floor(n);
123// n/=10;
124// }
125 // case where n approximates zero, set n to "positive" zero
126 if (abs(intp)==0)
127 {
128 if(abs(fract)<=RND_EPS)
129 {
130 n=0.000;
131 }
132 }
133 return n;
134}
135
136
137int main(int argc, char *argv[]) {
138
139 MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
140
141 try {
142
143 moab::Core mb_instance;
144 moab::Interface& moab = mb_instance;
145 int rank;
146 MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
147
148 //Reade parameters from line command
149 PetscBool flg = PETSC_TRUE;
150 char mesh_file_name[255];
151 ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
152 if(flg != PETSC_TRUE) {
153 SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
154 }
155 PetscInt order;
156 ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
157 if(flg != PETSC_TRUE) {
158 order = 5;
159 }
160
161 //Applied strain on the RVE (vector of length 6) strain=[xx, yy, zz, xy, xz, zy]^T
162 double myapplied_strain[6];
163 int nmax=6;
164 ierr = PetscOptionsGetRealArray(PETSC_NULL,"-myapplied_strain",myapplied_strain,&nmax,&flg); CHKERRQ(ierr);
165 VectorDouble applied_strain;
166 applied_strain.resize(6);
167 cblas_dcopy(6, &myapplied_strain[0], 1, &applied_strain(0), 1);
168 // cout<<"applied_strain ="<<applied_strain<<endl;
169
170
171 //Read mesh to MOAB
172 const char *option;
173 option = "";//"PARALLEL=BCAST;";//;DEBUG_IO";
174 rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
175 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
176 if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
177
178 //We need that for code profiling
179 PetscLogDouble t1,t2;
180 PetscLogDouble v1,v2;
181 ierr = PetscTime(&v1); CHKERRQ(ierr);
182 ierr = PetscGetCPUTime(&t1); CHKERRQ(ierr);
183
184 //Create MoFEM (Joseph) database
185 MoFEM::Core core(moab);
186 MoFEM::Interface& m_field = core;
187
188 //ref meshset ref level 0
189 ierr = m_field.seed_ref_level_3D(0,0); CHKERRQ(ierr);
190
191 // stl::bitset see for more details
192 BitRefLevel bit_level0;
193 bit_level0.set(0);
194 EntityHandle meshset_level0;
195 rval = moab.create_meshset(MESHSET_SET,meshset_level0); CHKERRQ_MOAB(rval);
196 ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
197 ierr = m_field.get_entities_by_ref_level(bit_level0,BitRefLevel().set(),meshset_level0); CHKERRQ(ierr);
198
199 //=======================================================================================================
200 //Add Periodic Prisims Between Triangles on -ve and +ve faces to implement periodic bounary conditions
201 //=======================================================================================================
202
203 //Populating the Multi-index container with -ve triangles
204 Range SurTrisNeg;
205 ierr = m_field.get_cubit_msId_entities_by_dimension(101,SIDESET,2,SurTrisNeg,true); CHKERRQ(ierr);
206 ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
207 Face_CenPos_Handle_multiIndex Face_CenPos_Handle_varNeg, Face_CenPos_Handle_varPos;
208 double TriCen[3], coords_Tri[9];
209 for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
210 const EntityHandle* conn_face; int num_nodes_Tri;
211
212 //get nodes attached to the face
213 rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
214 //get nodal coordinates
215 rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
216
217 //Find triangle centriod
218 TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
219 TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
220 TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
221
222 //round values to 3 disimal places
223 if(TriCen[0]>=0) TriCen[0]=double(int(TriCen[0]*1000.0+0.5))/1000.0; else TriCen[0]=double(int(TriCen[0]*1000.0-0.5))/1000.0; //-ve and +ve value
224 if(TriCen[1]>=0) TriCen[1]=double(int(TriCen[1]*1000.0+0.5))/1000.0; else TriCen[1]=double(int(TriCen[1]*1000.0-0.5))/1000.0;
225 if(TriCen[2]>=0) TriCen[2]=double(int(TriCen[2]*1000.0+0.5))/1000.0; else TriCen[2]=double(int(TriCen[2]*1000.0-0.5))/1000.0;
226
227 // cout<<"\n\n\nTriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
228 //fill the multi-index container with centriod coordinates and triangle handles
229 Face_CenPos_Handle_varNeg.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
230 // for(int ii=0; ii<3; ii++) cout<<"TriCen "<<TriCen[ii]<<endl;
231 // cout<<endl<<endl;
232 }
233
234 // double aaa;
235 // aaa=0.5011;
236 // cout<<"\n\n\n\n\nfloor(aaa+0.5) = "<<double(int(aaa*1000.0+0.5))/1000.0<<endl<<endl<<endl<<endl;
237 // aaa=-0.5011;
238 // cout<<"\n\n\n\n\nfloor(aaa+0.5) = "<<double(int(aaa*1000.0-0.5))/1000.0<<endl<<endl<<endl<<endl;
239
240
241 //Populating the Multi-index container with +ve triangles
242 Range SurTrisPos;
243 ierr = m_field.get_cubit_msId_entities_by_dimension(102,SIDESET,2,SurTrisPos,true); CHKERRQ(ierr);
244 ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 102 = %d\n",SurTrisPos.size()); CHKERRQ(ierr);
245 for(Range::iterator it = SurTrisPos.begin(); it!=SurTrisPos.end(); it++) {
246 const EntityHandle* conn_face; int num_nodes_Tri;
247
248 //get nodes attached to the face
249 rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
250 //get nodal coordinates
251 rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
252
253 //Find triangle centriod
254 TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
255 TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
256 TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
257
258 //round values to 3 disimal places
259 if(TriCen[0]>=0) TriCen[0]=double(int(TriCen[0]*1000.0+0.5))/1000.0; else TriCen[0]=double(int(TriCen[0]*1000.0-0.5))/1000.0;
260 if(TriCen[1]>=0) TriCen[1]=double(int(TriCen[1]*1000.0+0.5))/1000.0; else TriCen[1]=double(int(TriCen[1]*1000.0-0.5))/1000.0;
261 if(TriCen[2]>=0) TriCen[2]=double(int(TriCen[2]*1000.0+0.5))/1000.0; else TriCen[2]=double(int(TriCen[2]*1000.0-0.5))/1000.0;
262 // cout<<"\n\n\nTriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
263
264 //fill the multi-index container with centriod coordinates and triangle handles
265 Face_CenPos_Handle_varPos.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
266 }
267
268 //Find minimum and maximum X, Y and Z coordinates of the RVE (using multi-index container)
269 double XcoordMin, YcoordMin, ZcoordMin, XcoordMax, YcoordMax, ZcoordMax;
270 typedef Face_CenPos_Handle_multiIndex::index<xcoord_tag>::type::iterator Tri_Xcoord_iterator;
271 typedef Face_CenPos_Handle_multiIndex::index<ycoord_tag>::type::iterator Tri_Ycoord_iterator;
272 typedef Face_CenPos_Handle_multiIndex::index<zcoord_tag>::type::iterator Tri_Zcoord_iterator;
273 Tri_Xcoord_iterator XcoordMin_it, XcoordMax_it;
274 Tri_Ycoord_iterator YcoordMin_it, YcoordMax_it;
275 Tri_Zcoord_iterator ZcoordMin_it, ZcoordMax_it;
276
277 //XcoordMax_it-- because .end() will point iterator after the data range but .begin() will point the iteratore to the first value of range
278 XcoordMin_it=Face_CenPos_Handle_varNeg.get<xcoord_tag>().begin(); XcoordMin=XcoordMin_it->xcoord;
279 XcoordMax_it=Face_CenPos_Handle_varPos.get<xcoord_tag>().end(); XcoordMax_it--; XcoordMax=XcoordMax_it->xcoord;
280 YcoordMin_it=Face_CenPos_Handle_varNeg.get<ycoord_tag>().begin(); YcoordMin=YcoordMin_it->ycoord;
281 YcoordMax_it=Face_CenPos_Handle_varPos.get<ycoord_tag>().end(); YcoordMax_it--; YcoordMax=YcoordMax_it->ycoord;
282 ZcoordMin_it=Face_CenPos_Handle_varNeg.get<zcoord_tag>().begin(); ZcoordMin=ZcoordMin_it->zcoord;
283 ZcoordMax_it=Face_CenPos_Handle_varPos.get<zcoord_tag>().end(); ZcoordMax_it--; ZcoordMax=ZcoordMax_it->zcoord;
284
285// cout<<"XcoordMin "<<XcoordMin << " XcoordMax "<<XcoordMax <<endl;
286// cout<<"YcoordMin "<<YcoordMin << " YcoordMax "<<YcoordMax <<endl;
287// cout<<"ZcoordMin "<<ZcoordMin << " ZcoordMax "<<ZcoordMax <<endl;
288 double LxRVE, LyRVE, LzRVE;
289 LxRVE=XcoordMax-XcoordMin;
290 LyRVE=YcoordMax-YcoordMin;
291 LzRVE=ZcoordMax-ZcoordMin;
292
293
294 //Creating Prisims between triangles on -ve and +ve faces
295 typedef Face_CenPos_Handle_multiIndex::index<Tri_Hand_tag>::type::iterator Tri_Hand_iterator;
296 Tri_Hand_iterator Tri_Neg;
297 typedef Face_CenPos_Handle_multiIndex::index<Composite_xyzcoord>::type::iterator xyzcoord_iterator;
298 xyzcoord_iterator Tri_Pos;
299 Range PrismRange;
300 double XPos, YPos, ZPos;
301 //int count=0;
302
303 //loop over -ve triangles to create prisims elemenet between +ve and -ve triangles
304 for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
305
306 Tri_Neg=Face_CenPos_Handle_varNeg.get<Tri_Hand_tag>().find(*it);
307 //cout<<"xcoord= "<<Tri_iit->xcoord << " ycoord= "<< Tri_iit->ycoord << " ycoord= "<< Tri_iit->zcoord <<endl;
308
309 //corresponding +ve triangle
310 if(Tri_Neg->xcoord==XcoordMin){XPos=XcoordMax; YPos=Tri_Neg->ycoord; ZPos=Tri_Neg->zcoord;};
311 if(Tri_Neg->ycoord==YcoordMin){XPos=YPos=Tri_Neg->xcoord; YPos=YcoordMax; ZPos=Tri_Neg->zcoord;};
312 if(Tri_Neg->zcoord==ZcoordMin){XPos=YPos=Tri_Neg->xcoord; YPos=Tri_Neg->ycoord; ZPos=ZcoordMax; };
313 Tri_Pos=Face_CenPos_Handle_varPos.get<Composite_xyzcoord>().find(boost::make_tuple(XPos, YPos, ZPos));
314 // cout<<"Tri_Neg->xcoord= "<<Tri_Neg->xcoord << " Tri_Neg->ycoord "<< Tri_Neg->ycoord << " Tri_Neg->zcoord= "<< Tri_Neg->zcoord <<endl;
315 // cout<<"Tri_Pos->xcoord= "<<Tri_Pos->xcoord << " Tri_Pos->ycoord "<< Tri_Pos->ycoord << " Tri_Pos->zcoord= "<< Tri_Pos->zcoord <<endl;
316
317
318 //+ve and -ve nodes and their coords (+ve and -ve tiangles nodes can have matching problems, which can produce twisted prism)
319 EntityHandle PrismNodes[6];
320 vector<EntityHandle> TriNodesNeg, TriNodesPos;
321 double CoordNodeNeg[9], CoordNodePos[9];
322 rval = moab.get_connectivity(&(Tri_Neg->Tri_Hand),1,TriNodesNeg,true); CHKERRQ_MOAB(rval);
323 rval = moab.get_connectivity(&(Tri_Pos->Tri_Hand),1,TriNodesPos,true); CHKERRQ_MOAB(rval);
324 rval = moab.get_coords(&TriNodesNeg[0],3,CoordNodeNeg); MOAB_THROW(rval);
325 rval = moab.get_coords(&TriNodesPos[0],3,CoordNodePos); MOAB_THROW(rval);
326 for(int ii=0; ii<3; ii++){
327 PrismNodes[ii]=TriNodesNeg[ii];
328 }
329 // for(int ii=0; ii<3; ii++){
330 // cout<<"xcoord= "<<CoordNodeNeg[3*ii] << " ycoord= "<< CoordNodeNeg[3*ii+1] << " zcoord= "<< CoordNodeNeg[3*ii+2] <<endl;
331 // }
332 // for(int ii=0; ii<3; ii++){
333 // cout<<"xcoord= "<<CoordNodePos[3*ii] << " ycoord= "<< CoordNodePos[3*ii+1] << " zcoord= "<< CoordNodePos[3*ii+2] <<endl;
334 // }
335
336 //Match exact nodes to each other to avoide the problem of twisted prisms
337 double XNodeNeg, YNodeNeg, ZNodeNeg, XNodePos, YNodePos, ZNodePos;
338 for(int ii=0; ii<3; ii++){
339 if(Tri_Neg->xcoord==XcoordMin){XNodeNeg=XcoordMax; YNodeNeg=CoordNodeNeg[3*ii+1]; ZNodeNeg=CoordNodeNeg[3*ii+2];};
340 if(Tri_Neg->ycoord==YcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=YcoordMax; ZNodeNeg=CoordNodeNeg[3*ii+2];};
341 if(Tri_Neg->zcoord==ZcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=CoordNodeNeg[3*ii+1]; ZNodeNeg=ZcoordMax;};
342 for(int jj=0; jj<3; jj++){
343 //Round nodal coordinates to 3 dicimal places only for comparison
344 //round values to 3 disimal places
345 if(XNodeNeg>=0) XNodeNeg=double(int(XNodeNeg*1000.0+0.5))/1000.0; else XNodeNeg=double(int(XNodeNeg*1000.0-0.5))/1000.0;
346 if(YNodeNeg>=0) YNodeNeg=double(int(YNodeNeg*1000.0+0.5))/1000.0; else YNodeNeg=double(int(YNodeNeg*1000.0-0.5))/1000.0;
347 if(ZNodeNeg>=0) ZNodeNeg=double(int(ZNodeNeg*1000.0+0.5))/1000.0; else ZNodeNeg=double(int(ZNodeNeg*1000.0-0.5))/1000.0;
348
349 XNodePos=CoordNodePos[3*jj]; YNodePos=CoordNodePos[3*jj+1]; ZNodePos=CoordNodePos[3*jj+2];
350 if(XNodePos>=0) XNodePos=double(int(XNodePos*1000.0+0.5))/1000.0; else XNodePos=double(int(XNodePos*1000.0-0.5))/1000.0;
351 if(YNodePos>=0) YNodePos=double(int(YNodePos*1000.0+0.5))/1000.0; else YNodePos=double(int(YNodePos*1000.0-0.5))/1000.0;
352 if(ZNodePos>=0) ZNodePos=double(int(ZNodePos*1000.0+0.5))/1000.0; else ZNodePos=double(int(ZNodePos*1000.0-0.5))/1000.0;
353
354 if(XNodeNeg==XNodePos && YNodeNeg==YNodePos && ZNodeNeg==ZNodePos){
355 PrismNodes[3+ii]=TriNodesPos[jj];
356 break;
357 }
358 }
359 }
360 //prism nodes and their coordinates
361 double CoordNodesPrisms[18];
362 rval = moab.get_coords(PrismNodes,6,CoordNodesPrisms); MOAB_THROW(rval);
363 // for(int ii=0; ii<6; ii++){
364 // cout<<"xcoord= "<<CoordNodesPrisms[3*ii] << " ycoord= "<< CoordNodesPrisms[3*ii+1] << " zcoord= "<< CoordNodesPrisms[3*ii+2] <<endl;
365 // }
366 // cout<<endl<<endl;
367 //insertion of individula prism element and its addition to range PrismRange
368 EntityHandle PeriodicPrism;
369 rval = moab.create_element(MBPRISM,PrismNodes,6,PeriodicPrism); CHKERRQ_MOAB(rval);
370 PrismRange.insert(PeriodicPrism);
371
372 //Adding Prisims to Element LAGRANGE_ELEM (to loop over these prisims)
373 EntityHandle PrismRangeMeshset;
374 rval = moab.create_meshset(MESHSET_SET,PrismRangeMeshset); CHKERRQ_MOAB(rval);
375 rval = moab.add_entities(PrismRangeMeshset,PrismRange); CHKERRQ_MOAB(rval);
376 ierr = m_field.seed_ref_level_3D(PrismRangeMeshset,bit_level0); CHKERRQ(ierr);
377 ierr = m_field.get_entities_by_ref_level(bit_level0,BitRefLevel().set(),meshset_level0); CHKERRQ(ierr);
378
379 // //to see individual prisms
380 // Range Prism1;
381 // Prism1.insert(PeriodicPrism);
382 // EntityHandle out_meshset1;
383 // rval = moab.create_meshset(MESHSET_SET,out_meshset1); CHKERRQ_MOAB(rval);
384 // rval = moab.add_entities(out_meshset1,Prism1); CHKERRQ_MOAB(rval);
385 // ostringstream sss;
386 // sss << "Prism" << count << ".vtk"; count++;
387 // rval = moab.write_file(sss.str().c_str(),"VTK","",&out_meshset1,1); CHKERRQ_MOAB(rval);
388
389 }
390
391
392 //cout<<"PrismRange "<<PrismRange<<endl;
393 // //Saving prisms in interface.vtk
394 // EntityHandle out_meshset1;
395 // rval = moab.create_meshset(MESHSET_SET,out_meshset1); CHKERRQ_MOAB(rval);
396 // rval = moab.add_entities(out_meshset1,PrismRange); CHKERRQ_MOAB(rval);
397 // rval = moab.write_file("Prisms.vtk","VTK","",&out_meshset1,1); CHKERRQ_MOAB(rval);
398
399
400// //=======================================================================================================
401// //=======================================================================================================
402
403 //*
404 //Define problem
405
406 //Fields
407 int field_rank=3;
408 ierr = m_field.add_field("DISPLACEMENT",H1,AINSWORTH_LEGENDRE_BASE,field_rank); CHKERRQ(ierr);
409 ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,field_rank,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
410 ierr = m_field.add_field("LAGRANGE_MUL_PERIODIC",H1,AINSWORTH_LEGENDRE_BASE,field_rank); CHKERRQ(ierr); //For lagrange multipliers to control the periodic motion
411 ierr = m_field.add_field("LAGRANGE_MUL_RIGID_TRANS",NOFIELD,NOBASE,3); CHKERRQ(ierr); //To control the rigid body motion (3 Traslations and 3 rotations)
412
413 //add entitities (by tets) to the field
414 //============================================================================================================
415 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"DISPLACEMENT"); CHKERRQ(ierr);
416 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
417
418 //Adding only -ve surfaces to the field LAGRANGE_MUL_PERIODIC (in periodic boundary conditions size of C (3M/2 x 3N))
419 //to create meshset from range SurTrisNeg
420 EntityHandle SurTrisNegMeshset;
421 rval = moab.create_meshset(MESHSET_SET,SurTrisNegMeshset); CHKERRQ_MOAB(rval);
422 rval = moab.add_entities(SurTrisNegMeshset,SurTrisNeg); CHKERRQ_MOAB(rval);
423 ierr = m_field.add_ents_to_field_by_type(SurTrisNegMeshset,MBTRI,"LAGRANGE_MUL_PERIODIC"); CHKERRQ(ierr);
424 //============================================================================================================
425
426 ierr = m_field.set_field_order(0,MBTET,"DISPLACEMENT",order); CHKERRQ(ierr);
427 ierr = m_field.set_field_order(0,MBTRI,"DISPLACEMENT",order); CHKERRQ(ierr);
428 ierr = m_field.set_field_order(0,MBEDGE,"DISPLACEMENT",order); CHKERRQ(ierr);
429 ierr = m_field.set_field_order(0,MBVERTEX,"DISPLACEMENT",1); CHKERRQ(ierr);
430
431 ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
432 ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
433 ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
434 ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
435
436 //ierr = m_field.set_field_order(0,MBPRISM,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
437 ierr = m_field.set_field_order(0,MBTRI,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
438 ierr = m_field.set_field_order(0,MBEDGE,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
439 ierr = m_field.set_field_order(0,MBVERTEX,"LAGRANGE_MUL_PERIODIC",1); CHKERRQ(ierr);
440
441 //Define FE
442 //define eleatic element
443 boost::shared_ptr<Hooke<adouble> > hooke_adouble(new Hooke<adouble>);
444 boost::shared_ptr<Hooke<double> > hooke_double(new Hooke<double>);
445 NonlinearElasticElement elastic(m_field,2);
446 ierr = elastic.setBlocks(hooke_double,hooke_adouble); CHKERRQ(ierr);
447 ierr = elastic.addElement("ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
448 ierr = elastic.setOperators("DISPLACEMENT","MESH_NODE_POSITIONS",false,true); CHKERRQ(ierr);
449
450 BCs_RVELagrange_Periodic lagrangian_element_periodic(m_field);
451 BCs_RVELagrange_Trac lagrangian_element_trac(m_field);
452
453 lagrangian_element_periodic.addLagrangiangElement("LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",PrismRange);
454 lagrangian_element_trac.addLagrangiangElement("LAGRANGE_ELEM_TRANS","DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS","MESH_NODE_POSITIONS");
455
456 //define problems
457 ierr = m_field.add_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
458
459 //set finite elements for problem
460 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","ELASTIC"); CHKERRQ(ierr);
461 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM"); CHKERRQ(ierr);
462 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
463
464
465 //set refinement level for problem
466 ierr = m_field.modify_problem_ref_level_add_bit("ELASTIC_MECHANICS",bit_level0); CHKERRQ(ierr);
467
468 /***/
469 //Declare problem
470 /****/
471 //build database
472
473 //build field
474 ierr = m_field.build_fields(); CHKERRQ(ierr);
475 Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
476 ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
477
478
479 //build finite elemnts
480 ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
481
482 //build adjacencies
483 ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
484
485 //build problem
486 ierr = m_field.build_problems(); CHKERRQ(ierr);
487
488 /****/
489 //mesh partitioning
490
491 //partition
492 ierr = m_field.partition_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
493 ierr = m_field.partition_finite_elements("ELASTIC_MECHANICS"); CHKERRQ(ierr);
494 //what are ghost nodes, see Petsc Manual
495 ierr = m_field.partition_ghost_dofs("ELASTIC_MECHANICS"); CHKERRQ(ierr);
496
497 //print bcs
498 ierr = m_field.print_cubit_displacement_set(); CHKERRQ(ierr);
499 ierr = m_field.print_cubit_force_set(); CHKERRQ(ierr);
500 //print block sets with materials
501 ierr = m_field.print_cubit_materials_set(); CHKERRQ(ierr);
502
503 //create matrices (here F, D and Aij are matrices for the full problem)
504 Vec D;
505 vector<Vec> F(6);
506 ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",ROW,&F[0]); CHKERRQ(ierr);
507 for(int ii = 1;ii<6;ii++) {
508 ierr = VecDuplicate(F[0],&F[ii]); CHKERRQ(ierr);
509 }
510
511 ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",COL,&D); CHKERRQ(ierr);
512 Mat Aij;
513 ierr = m_field.MatCreateMPIAIJWithArrays("ELASTIC_MECHANICS",&Aij); CHKERRQ(ierr);
514
515 for(int ii = 0;ii<6;ii++) {
516 ierr = VecZeroEntries(F[ii]); CHKERRQ(ierr);
517 ierr = VecGhostUpdateBegin(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
518 ierr = VecGhostUpdateEnd(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
519 }
520
521 ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
522
523
524 //Assemble F and Aij
525 elastic.getLoopFeLhs().snes_B = Aij;
526 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC",elastic.getLoopFeLhs()); CHKERRQ(ierr);
527
528 lagrangian_element_periodic.setRVEBCsOperators("DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",Aij,F);
529 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
530 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCRhs()); CHKERRQ(ierr);
531
532 BCs_RVELagrange_Trac_Rigid_Trans lagrangian_element_rigid_body_trans(m_field);
533 lagrangian_element_rigid_body_trans.setRVEBCsRigidBodyTranOperators("DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS",Aij, lagrangian_element_periodic.setOfRVEBC);
534 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS",lagrangian_element_rigid_body_trans.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
535
536 for(int ii = 0;ii<6;ii++) {
537 ierr = VecGhostUpdateBegin(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
538 ierr = VecGhostUpdateEnd(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
539 ierr = VecAssemblyBegin(F[ii]); CHKERRQ(ierr);
540 ierr = VecAssemblyEnd(F[ii]); CHKERRQ(ierr);
541 }
542
543 ierr = MatAssemblyBegin(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
544 ierr = MatAssemblyEnd(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
545
546 // Volume calculation
547 //==============================================================================================================================
548 VolumeElementForcesAndSourcesCore vol_elem(m_field);
549 Vec volume_vec;
550 int volume_vec_ghost[] = { 0 };
551 ierr = VecCreateGhost(
552 PETSC_COMM_WORLD,(!m_field.get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
553 ); CHKERRQ(ierr);
554 ierr = VecZeroEntries(volume_vec); CHKERRQ(ierr);
555 vol_elem.getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
556
557 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC", vol_elem); CHKERRQ(ierr);
558
559 ierr = VecAssemblyBegin(volume_vec); CHKERRQ(ierr);
560 ierr = VecAssemblyEnd(volume_vec); CHKERRQ(ierr);
561 double rve_volume;
562 ierr = VecSum(volume_vec,&rve_volume); CHKERRQ(ierr);
563
564 ierr = VecView(volume_vec,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
565 ierr = PetscPrintf(PETSC_COMM_WORLD,"RVE Volume %3.2g\n",rve_volume); CHKERRQ(ierr);
566 ierr = VecDestroy(&volume_vec);
567
568
569 //==============================================================================================================================
570 //Post processing
571 //==============================================================================================================================
573 bool doPreProcess;
574 bool doPostProcess;
577 doPreProcess(true),
578 doPostProcess(true)
579 {}
580 void setDoPreProcess() { doPreProcess = true; }
581 void unSetDoPreProcess() { doPreProcess = false; }
582 void setDoPostProcess() { doPostProcess = true; }
583 void unSetDoPostProcess() { doPostProcess = false; }
584
585
586
587 PetscErrorCode preProcess() {
588 PetscFunctionBegin;
589 if(doPreProcess) {
591 }
592 PetscFunctionReturn(0);
593 }
594 PetscErrorCode postProcess() {
595 PetscFunctionBegin;
596 if(doPostProcess) {
598 }
599 PetscFunctionReturn(0);
600 }
601 };
602
603 MyPostProc post_proc(m_field);
604
605 ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
606 ierr = post_proc.addFieldValuesPostProc("DISPLACEMENT"); CHKERRQ(ierr);
607 ierr = post_proc.addFieldValuesGradientPostProc("DISPLACEMENT"); CHKERRQ(ierr);
608 ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
609 ierr = post_proc.addFieldValuesGradientPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
610
611 for(
612 map<int,NonlinearElasticElement::BlockData>::iterator sit = elastic.setOfBlocks.begin();
613 sit != elastic.setOfBlocks.end(); sit++
614 ) {
615 post_proc.getOpPtrVector().push_back(
616 new PostPorcStress(
617 post_proc.postProcMesh,
618 post_proc.mapGaussPts,
619 "DISPLACEMENT",
620 sit->second,
621 post_proc.commonData,
622 false
623 )
624 );
625 }
626 ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
627 ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
628 ierr = m_field.set_global_ghost_vector(
629 "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
630 ); CHKERRQ(ierr);
631
632 //==============================================================================================================================
633
634 //Solver
635 KSP solver;
636 ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
637 ierr = KSPSetOperators(solver,Aij,Aij); CHKERRQ(ierr);
638 ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
639 ierr = KSPSetUp(solver); CHKERRQ(ierr);
640
641 MatrixDouble Dmat;
642 Dmat.resize(6,6);
643 Dmat.clear();
644
645 //calculate honmogenised stress
646 //create a vector for 6 components of homogenized stress
647 Vec stress_homo;
648 int stress_homo_ghost[] = { 0,1,2,3,4,5,6 };
649 ierr = VecCreateGhost(
650 PETSC_COMM_WORLD,(!m_field.get_comm_rank())?6:0,6,6,stress_homo_ghost,&stress_homo
651 ); CHKERRQ(ierr);
652
653 lagrangian_element_periodic.setRVEBCsHomoStressOperators("DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",stress_homo);
654
655 PetscScalar *avec;
656 ierr = VecGetArray(stress_homo,&avec); CHKERRQ(ierr);
657 for(int ii = 0;ii<6;ii++) {
658 ierr = VecZeroEntries(D); CHKERRQ(ierr);
659 ierr = KSPSolve(solver,F[ii],D); CHKERRQ(ierr);
660 ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
661 ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
662 ierr = m_field.set_global_ghost_vector(
663 "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
664 ); CHKERRQ(ierr);
665
666 //post-processing the resutls
667 post_proc.setDoPreProcess();
668 post_proc.unSetDoPostProcess();
669 ierr = m_field.loop_finite_elements(
670 "ELASTIC_MECHANICS","ELASTIC",post_proc
671 ); CHKERRQ(ierr);
672 post_proc.unSetDoPreProcess();
673 post_proc.setDoPostProcess();
674 {
675 ostringstream sss;
676 sss << "mode_" << "Disp" << "_" << ii << ".h5m";
677 rval = post_proc.postProcMesh.write_file(
678 sss.str().c_str(),"MOAB","PARALLEL=WRITE_PART"
679 ); CHKERRQ_MOAB(rval);
680 }
681
682 ierr = VecZeroEntries(stress_homo); CHKERRQ(ierr);
683 ierr = m_field.loop_finite_elements(
684 "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCStress()
685 ); CHKERRQ(ierr);
687 PETSC_NULL,"-my_rve_volume",&rve_volume,PETSC_NULL
688 ); CHKERRQ(ierr);
689 ierr = VecAssemblyBegin(stress_homo); CHKERRQ(ierr);
690 ierr = VecAssemblyEnd(stress_homo); CHKERRQ(ierr);
691 ierr = VecScale(stress_homo,1.0/rve_volume); CHKERRQ(ierr);
692 ierr = VecGhostUpdateBegin(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
693 ierr = VecGhostUpdateEnd(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
694 for(int jj=0; jj<6; jj++) {
695 Dmat(jj,ii) = avec[jj];
696 }
697 }
698
699 ierr = VecRestoreArray(stress_homo,&avec); CHKERRQ(ierr);
700 PetscPrintf(
701 PETSC_COMM_WORLD,"\nHomogenised Stiffens Matrix = \n\n"
702 );
703
704 for(int ii=0; ii<6; ii++) {
705 PetscPrintf(
706 PETSC_COMM_WORLD,
707 "stress %d\t\t%8.5e\t\t%8.5e\t\t%8.5e\t\t%8.5e\t\t%8.5e\t\t%8.5e\n",
708 ii,Dmat(ii,0),Dmat(ii,1),Dmat(ii,2),Dmat(ii,3),Dmat(ii,4),Dmat(ii,5)
709 );
710 }
711 //==============================================================================================================================
712
713 //Open mesh_file_name.txt for writing
714 ofstream myfile;
715 myfile.open ((string(mesh_file_name)+".txt").c_str());
716
717 //Output displacements
718 myfile << "<<<< Homonenised stress >>>>>" << endl;
719
720 if(pcomm->rank()==0){
721 for(int ii=0; ii<6; ii++){
722 myfile << boost::format("%.3lf") % roundn(Dmat(ii,0)) << endl;
723 avec++;
724 }
725 }
726 //Close mesh_file_name.txt
727 myfile.close();
728
729
730
731 //detroy matrices
732 for(int ii = 0;ii<6;ii++) {
733 ierr = VecDestroy(&F[ii]); CHKERRQ(ierr);
734 }
735 ierr = VecDestroy(&D); CHKERRQ(ierr);
736 ierr = MatDestroy(&Aij); CHKERRQ(ierr);
737 ierr = KSPDestroy(&solver); CHKERRQ(ierr);
738 ierr = VecDestroy(&stress_homo); CHKERRQ(ierr);
739
740 ierr = PetscTime(&v2);CHKERRQ(ierr);
741 ierr = PetscGetCPUTime(&t2);CHKERRQ(ierr);
742
743 PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Time = %f CPU Time = %f\n",pcomm->rank(),v2-v1,t2-t1);
744
745 }
747
749
750}
Implementation of linear elastic material.
Operators and data structures for non-linear elastic analysis.
Post-process fields on refined mesh.
Post-processing stresses for non-linear analysis.
DEPRECATED typedef PostProcStress PostPorcStress
Operator can be used with any volume element to calculate sum of volumes of all volumes in the set.
@ COL
Definition: definitions.h:123
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ MF_ZERO
Definition: definitions.h:98
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ NOBASE
Definition: definitions.h:59
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:541
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
@ H1
continuous field
Definition: definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
@ SIDESET
Definition: definitions.h:147
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:454
FTensor::Index< 'n', SPACE_DIM > n
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string &name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
multi_index_container< Face_CenPos_Handle, indexed_by< ordered_non_unique< tag< xcoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::xcoord > >, ordered_non_unique< tag< ycoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::ycoord > >, ordered_non_unique< tag< zcoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::zcoord > >, ordered_unique< tag< Tri_Hand_tag >, member< Face_CenPos_Handle, const EntityHandle,&Face_CenPos_Handle::Tri_Hand > >, ordered_unique< tag< Composite_xyzcoord >, composite_key< Face_CenPos_Handle, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::xcoord >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::ycoord >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::zcoord > > > > > Face_CenPos_Handle_multiIndex
int main(int argc, char *argv[])
double roundn(double n)
static char help[]
char mesh_file_name[255]
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
CoreTmp< 0 > Core
Definition: Core.hpp:1086
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
const double D
diffusivity
map< int, RVEBC_Data > setOfRVEBC
maps side set id with appropriate FluxData
PetscErrorCode addLagrangiangElement(const string element_name, const string field_name, const string lagrang_field_name, const string mesh_nodals_positions)
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec stress_homo)
PetscErrorCode addLagrangiangElement(const string element_name, const string field_name, const string lagrang_field_name, const string mesh_nodals_positions, Range &periodic_prisms)
PetscErrorCode setRVEBCsOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &f)
PetscErrorCode setRVEBCsRigidBodyTranOperators(string field_name, string lagrang_field_name, Mat aij, map< int, RVEBC_Data > setOfRVEBC)
Face_CenPos_Handle(double _xcoord, double _ycoord, double _zcoord, const EntityHandle _Tri_Hand)
Hook equation.
Definition: Hooke.hpp:21
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual int get_comm_rank() const =0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
DEPRECATED MoFEMErrorCode seed_ref_level_3D(const EntityHandle meshset, const BitRefLevel &bit, int verb=-1)
seed 2D entities in the meshset and their adjacencies (only TETs adjacencies) in a particular BitRefL...
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Projection of edge entities with one mid-node on hierarchical basis.
Mat & snes_B
preconditioner of jacobian matrix
Definition: plot_base.cpp:34
MoFEMErrorCode postProcess()
Definition: plot_base.cpp:622
MoFEMErrorCode generateReferenceElementMesh()
Definition: plot_base.cpp:397
MoFEMErrorCode preProcess()
Definition: plot_base.cpp:612
structure grouping operators and data used for calculation of nonlinear elastic element
MyVolumeFE & getLoopFeLhs()
get lhs volume element
MoFEMErrorCode addElement(const std::string element_name, const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false)
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
MoFEMErrorCode setBlocks(boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double > > materialDoublePtr, boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble > > materialAdoublePtr)
MoFEMErrorCode setOperators(const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false, const bool field_disp=false)
Set operators to calculate left hand tangent matrix and right hand residual.
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode postProcess()
function is run at the end of loop
Post processing.
Calculate volume.