v0.14.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 
19 #include <Projection10NodeCoordsOnField.hpp>
20 #include <petsctime.h>
21 
22 // #include <FEMethod_LowLevelStudent.hpp>
23 // #include <FEMethod_UpLevelStudent.hpp>
24 
25 #include <adolc/adolc.h>
26 #include <NonLinearElasticElement.hpp>
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 
36 #include <MethodForForceScaling.hpp>
37 #include <TimeForceScale.hpp>
38 #include <VolumeCalculation.hpp>
39 #include <BCs_RVELagrange_Disp.hpp>
40 #include <BCs_RVELagrange_Trac.hpp>
41 #include "BCs_RVELagrange_Periodic.hpp"
42 #include <BCs_RVELagrange_Trac_Rigid_Trans.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 
61 struct Face_CenPos_Handle {
62  double xcoord, ycoord, zcoord;
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 
68 struct xcoord_tag {}; //tags to used in the multindex container
69 struct ycoord_tag {};
70 struct zcoord_tag {};
71 struct Tri_Hand_tag {};
72 struct Composite_xyzcoord {};
73 
74 typedef multi_index_container<
75 Face_CenPos_Handle,
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<
92 Face_CenPos_Handle,
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> > >
96 > > Face_CenPos_Handle_multiIndex;
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):
576  PostProcVolumeOnRefinedMesh(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) {
590  ierr = PostProcVolumeOnRefinedMesh::preProcess(); CHKERRQ(ierr);
591  }
592  PetscFunctionReturn(0);
593  }
594  PetscErrorCode postProcess() {
595  PetscFunctionBegin;
596  if(doPostProcess) {
597  ierr = PostProcVolumeOnRefinedMesh::postProcess(); CHKERRQ(ierr);
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);
686  ierr = PetscOptionsGetReal(
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 
748  MoFEM::Core::Finalize();
749 
750 }
SIDESET
@ SIDESET
Definition: definitions.h:160
MoFEM::CoreInterface::loop_finite_elements
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.
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::CoreInterface::loop_dofs
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.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
RND_EPS
#define RND_EPS
Definition: homogenisation_periodic_atom_test.cpp:102
help
static char help[]
Definition: homogenisation_periodic_atom_test.cpp:54
BCs_RVELagrange_Periodic
Definition: BCs_RVELagrange_Periodic.hpp:19
MyPostProc
Definition: plot_base.cpp:34
BCs_RVELagrange_Disp.hpp
EntityHandle
Face_CenPos_Handle_multiIndex
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
Definition: homogenisation_periodic_atom_test.cpp:96
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:467
roundn
double roundn(double n)
Definition: homogenisation_periodic_atom_test.cpp:105
NOBASE
@ NOBASE
Definition: definitions.h:59
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
PostProcTemplateVolumeOnRefinedMesh< MoFEM::VolumeElementForcesAndSourcesCore >::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: PostProcOnRefMesh.hpp:797
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
Face_CenPos_Handle::Tri_Hand
const EntityHandle Tri_Hand
Definition: homogenisation_periodic_atom_test.cpp:63
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM.hpp
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MOAB_THROW
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:554
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::CoreInterface::add_ents_to_field_by_type
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.
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
BCs_RVELagrange_Periodic::setRVEBCsHomoStressOperators
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec stress_homo)
Definition: BCs_RVELagrange_Periodic.hpp:1312
Composite_xyzcoord
Definition: homogenisation_periodic_atom_test.cpp:72
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
xcoord_tag
Definition: homogenisation_periodic_atom_test.cpp:68
PostProcStresses.hpp
Post-processing stresses for non-linear analysis.
PostProcTemplateVolumeOnRefinedMesh< MoFEM::VolumeElementForcesAndSourcesCore >::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: PostProcOnRefMesh.hpp:696
main
int main(int argc, char *argv[])
Definition: homogenisation_periodic_atom_test.cpp:137
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
BCs_RVELagrange_Periodic::addLagrangiangElement
PetscErrorCode addLagrangiangElement(const string element_name, const string field_name, const string lagrang_field_name, const string mesh_nodals_positions, Range &periodic_prisms)
Definition: BCs_RVELagrange_Periodic.hpp:52
MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:152
BCs_RVELagrange_Disp::getLoopFeRVEBCLhs
MyTriFE & getLoopFeRVEBCLhs()
Definition: BCs_RVELagrange_Disp.hpp:35
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
Hooke.hpp
Implementation of linear elastic material.
double
BCs_RVELagrange_Trac_Rigid_Trans.hpp
COL
@ COL
Definition: definitions.h:136
NonLinearElasticElement.hpp
Operators and data structures for non-linear elastic analysis.
Projection10NodeCoordsOnField.hpp
BCs_RVELagrange_Trac
Definition: BCs_RVELagrange_Trac.hpp:22
PostProcOnRefMesh.hpp
Post-process fields on refined mesh.
NonlinearElasticElement::addElement
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)
Definition: NonLinearElasticElement.cpp:1120
VolumeCalculation.hpp
Operator can be used with any volume element to calculate sum of volumes of all volumes in the set.
Face_CenPos_Handle
Definition: homogenisation_periodic_atom_test.cpp:61
MyPostProc::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Definition: plot_base.cpp:427
BCs_RVELagrange_Periodic::setRVEBCsOperators
PetscErrorCode setRVEBCsOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &f)
Definition: BCs_RVELagrange_Periodic.hpp:1176
BCs_RVELagrange_Periodic::getLoopFeRVEBCRhs
MyPrismFE & getLoopFeRVEBCRhs()
Definition: BCs_RVELagrange_Periodic.hpp:39
BCs_RVELagrange_Trac_Rigid_Trans::setRVEBCsRigidBodyTranOperators
PetscErrorCode setRVEBCsRigidBodyTranOperators(string field_name, string lagrang_field_name, Mat aij, map< int, RVEBC_Data > setOfRVEBC)
Definition: BCs_RVELagrange_Trac_Rigid_Trans.hpp:109
BCs_RVELagrange_Periodic::getLoopFeRVEBCLhs
MyPrismFE & getLoopFeRVEBCLhs()
Definition: BCs_RVELagrange_Periodic.hpp:38
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
MoFEM::CoreInterface::modify_problem_ref_level_add_bit
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
Face_CenPos_Handle::Face_CenPos_Handle
Face_CenPos_Handle(double _xcoord, double _ycoord, double _zcoord, const EntityHandle _Tri_Hand)
Definition: homogenisation_periodic_atom_test.cpp:64
NonlinearElasticElement::setOperators
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.
Definition: NonLinearElasticElement.cpp:1153
BCs_RVELagrange_Disp::setOfRVEBC
map< int, RVEBC_Data > setOfRVEBC
maps side set id with appropriate FluxData
Definition: BCs_RVELagrange_Disp.hpp:56
BCs_RVELagrange_Trac_Rigid_Trans
Definition: BCs_RVELagrange_Trac_Rigid_Trans.hpp:18
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:22
Face_CenPos_Handle::zcoord
double zcoord
Definition: homogenisation_periodic_atom_test.cpp:62
convert.n
n
Definition: convert.py:82
Range
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition: DeprecatedPetsc.hpp:192
MoFEM::CoreTmp< 0 >::Initialize
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
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:955
PostPorcStress
DEPRECATED typedef PostProcStress PostPorcStress
Definition: PostProcStresses.hpp:193
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
BCs_RVELagrange_Periodic.hpp
NonlinearElasticElement::setOfBlocks
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
Definition: NonLinearElasticElement.hpp:100
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
Tri_Hand_tag
Definition: homogenisation_periodic_atom_test.cpp:71
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::DeprecatedCoreInterface::seed_ref_level_3D
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...
Definition: DeprecatedCoreInterface.cpp:18
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
BCs_RVELagrange_Disp::addLagrangiangElement
PetscErrorCode addLagrangiangElement(const string element_name, const string field_name, const string lagrang_field_name, const string mesh_nodals_positions)
Definition: BCs_RVELagrange_Disp.hpp:58
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::CoreInterface::modify_problem_add_finite_element
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
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::CoreInterface::set_field_order
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.
VolumeCalculation
Calculate volume.
Definition: VolumeCalculation.hpp:17
BCs_RVELagrange_Trac.hpp
MoFEM::CoreInterface::add_problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEM::CoreInterface::add_field
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.
ycoord_tag
Definition: homogenisation_periodic_atom_test.cpp:69
NonlinearElasticElement::setBlocks
MoFEMErrorCode setBlocks(boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double >> materialDoublePtr, boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble >> materialAdoublePtr)
Definition: NonLinearElasticElement.cpp:1086
BCs_RVELagrange_Periodic::getLoopFeRVEBCStress
MyPrismFE & getLoopFeRVEBCStress()
Definition: BCs_RVELagrange_Periodic.hpp:42
F
@ F
Definition: free_surface.cpp:394
Hooke
Hook equation.
Definition: Hooke.hpp:19
NOFIELD
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
zcoord_tag
Definition: homogenisation_periodic_atom_test.cpp:70
Generated by Doxygen 1.8.17 and hosted at University of Glasgow