v0.13.0
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>
17 using 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>
39 #include <BCs_RVELagrange_Disp.hpp>
40 #include <BCs_RVELagrange_Trac.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 
54 static 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 
62  double xcoord, ycoord, zcoord;
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 
68 struct xcoord_tag {}; //tags to used in the multindex container
69 struct ycoord_tag {};
70 struct zcoord_tag {};
71 struct Tri_Hand_tag {};
73 
74 typedef multi_index_container<
76 indexed_by<
77 ordered_non_unique<
78 tag<xcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord> >,
79 
80 ordered_non_unique<
81 tag<ycoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord> >,
82 
83 ordered_non_unique<
84 tag<zcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> >,
85 
86 ordered_unique<
87 tag<Tri_Hand_tag>, member<Face_CenPos_Handle,const EntityHandle,&Face_CenPos_Handle::Tri_Hand> >,
88 
89 ordered_unique<
90 tag<Composite_xyzcoord>,
91 composite_key<
93 member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord>,
94 member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord>,
95 member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> > >
97 
98 //============================================================================================
99 //============================================================================================
100 
101 
102 #define RND_EPS 1e-6
103 
104 //Rounding
105 double 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 
137 int 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  //==============================================================================================================================
572  struct MyPostProc: public PostProcVolumeOnRefinedMesh {
573  bool doPreProcess;
574  bool doPostProcess;
575  MyPostProc(MoFEM::Interface &m_field):
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  }
746  CATCH_ERRORS;
747 
749 
750 }
static Index< 'n', 3 > n
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:136
@ ROW
Definition: definitions.h:136
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ MF_ZERO
Definition: definitions.h:111
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ NOBASE
Definition: definitions.h:72
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:554
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:97
@ H1
continuous field
Definition: definitions.h:98
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
@ SIDESET
Definition: definitions.h:160
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:467
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:87
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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:1096
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
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:33
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:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
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:52
MoFEMErrorCode generateReferenceElementMesh()
Definition: plot_base.cpp:354
structure grouping operators and data used for calculation of nonlinear elastic 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)
MyVolumeFE & getLoopFeLhs()
get lhs volume element
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.