v0.14.0
rve_mechanical.cpp
Go to the documentation of this file.
1 /** \file rve_mechanical.cpp
2  * \brief Calculates stiffness matrix for elastic RVE.
3 
4  Three types of boundary conditions are implemented, i.e.
5  HOMOBCDISP, HOMOBCPERIODIC, HOMOBCTRAC, NITSCHE.
6 
7  NITSHCE method allow to apply periodic boundary conditions
8  to arbitrary convex shape RVE.
9 
10  */
11 
12 
13 
14 /* This file is part of MoFEM.
15  * MoFEM is free software: you can redistribute it and/or modify it under
16  * the terms of the GNU Lesser General Public License as published by the
17  * Free Software Foundation, either version 3 of the License, or (at your
18  * option) any later version.
19  *
20  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
21  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
22  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
23  * License for more details.
24  *
25  * You should have received a copy of the GNU Lesser General Public
26  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
27 
28 #include <MoFEM.hpp>
29 using namespace MoFEM;
30 
31 #include <Projection10NodeCoordsOnField.hpp>
32 #include <petsctime.h>
33 
34 #include <boost/numeric/ublas/vector_proxy.hpp>
35 
36 #include <MethodForForceScaling.hpp>
37 #include <TimeForceScale.hpp>
38 
39 #include <BCs_RVELagrange_Disp.hpp>
40 #include <BCs_RVELagrange_Trac.hpp>
41 #include <BCs_RVELagrange_Trac_Rigid_Rot.hpp>
42 #include <BCs_RVELagrange_Trac_Rigid_Trans.hpp>
43 #include <BCs_RVELagrange_Periodic.hpp>
44 
45 #ifndef WITH_ADOL_C
46  #error "MoFEM need to be compiled with ADOL-C"
47 #endif
48 
49 #include <adolc/adolc.h>
50 #include <NonLinearElasticElement.hpp>
51 #include <Hooke.hpp>
52 #include <boost/numeric/ublas/symmetric.hpp>
53 #include <SmallTransverselyIsotropic.hpp>
54 #include <VolumeCalculation.hpp>
55 
56 #include <boost/ptr_container/ptr_map.hpp>
57 #include <PostProcOnRefMesh.hpp>
58 #include <PostProcStresses.hpp>
59 
60 extern "C" {
61  void tetcircumcenter_tp(double a[3],double b[3],double c[3], double d[3],
62  double circumcenter[3],double *xi,double *eta,double *zeta);
63  void tricircumcenter3d_tp(double a[3],double b[3],double c[3],
64  double circumcenter[3],double *xi,double *eta);
65  //#include <triangle_ncc_rule.h>
66 }
67 
68 // #include <NitscheMethod.hpp>
69 #include <moab/AdaptiveKDTree.hpp>
70 // #include <NitschePeriodicMethod.hpp>
71 
72 #include <algorithm>
73 
74 
75 
76 
77 static char help[] = "...\n\n";
78 
79 enum HomoBCTypes {
80  HOMOBCDISP,
81  HOMOBCPERIODIC,
82  HOMOBCTRAC,
83  NITSCHE,
84  NOHOMOBC
85 };
86 
87 const char *homo_bc_names[] = {
88  "disp",
89  "periodic",
90  "trac",
91  "nitsche",
92  "nohomobc"
93 };
94 
95 //=================================================================================================================================
96 //Define class and multindex container to store data for traiangles on the boundary of the RVE (it cannot be defined within main)
97 //=================================================================================================================================
98 
99 struct Face_CenPos_Handle {
100  double xcoord, ycoord, zcoord;
101  const EntityHandle Tri_Hand;
102  Face_CenPos_Handle(double _xcoord, double _ycoord, double _zcoord, const EntityHandle _Tri_Hand):xcoord(_xcoord),
103  ycoord(_ycoord), zcoord(_zcoord), Tri_Hand(_Tri_Hand) {}
104 };
105 
106 struct xcoord_tag {}; //tags to used in the multindex container
107 struct ycoord_tag {};
108 struct zcoord_tag {};
109 struct Tri_Hand_tag {};
110 struct Composite_xyzcoord {};
111 
112 typedef multi_index_container<
113 Face_CenPos_Handle,
114 indexed_by<
115 ordered_non_unique<
116 tag<xcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord> >,
117 
118 ordered_non_unique<
119 tag<ycoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord> >,
120 
121 ordered_non_unique<
122 tag<zcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> >,
123 
124 ordered_unique<
125 tag<Tri_Hand_tag>, member<Face_CenPos_Handle,const EntityHandle,&Face_CenPos_Handle::Tri_Hand> >,
126 
127 ordered_unique<
128 tag<Composite_xyzcoord>,
129 composite_key<
130 Face_CenPos_Handle,
131 member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord>,
132 member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord>,
133 member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> > >
134 > > Face_CenPos_Handle_multiIndex;
135 
136 int main(int argc, char *argv[]) {
137 
138  MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
139 
140  try {
141 
142  moab::Core mb_instance;
143  moab::Interface& moab = mb_instance;
144  int rank;
145  MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
146 
147  //Reade parameters from line command
148  PetscBool flg = PETSC_TRUE;
149  char mesh_file_name[255];
150  ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
151  if(flg != PETSC_TRUE) {
152  SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_FOUND,"*** ERROR -my_file (MESH FILE NEEDED)");
153  }
154 
155  PetscInt order;
156  ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
157  if(flg != PETSC_TRUE) {
158  order = 1;
159  }
160 
161  PetscInt choise_value = NOHOMOBC;
162  ierr = PetscOptionsGetEList(
163  NULL,"-my_bc_type",homo_bc_names,NOHOMOBC,&choise_value,&flg
164  ); CHKERRQ(ierr);
165  if(flg != PETSC_TRUE) {
166  SETERRQ(PETSC_COMM_SELF,MOFEM_IMPOSIBLE_CASE,"*** Boundary conditions not set");
167  }
168 
169  //Read mesh to MOAB
170  const char *option;
171  option = "";//"PARALLEL=BCAST;";//;DEBUG_IO";
172  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
173  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
174  if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
175 
176  //We need that for code profiling
177  PetscLogDouble t1,t2;
178  PetscLogDouble v1,v2;
179  ierr = PetscTime(&v1); CHKERRQ(ierr);
180  ierr = PetscGetCPUTime(&t1); CHKERRQ(ierr);
181 
182  //Create MoFEM (Joseph) database
183  MoFEM::Core core(moab);
184  MoFEM::Interface& m_field = core;
185 
186  vector<BitRefLevel> bit_levels;
187  {
188  Tag th_meshset_info;
189  int def_meshset_info[2] = {0,0};
190  rval = moab.tag_get_handle(
191  "MESHSET_INFO",2,MB_TYPE_INTEGER,th_meshset_info,MB_TAG_CREAT|MB_TAG_SPARSE,&def_meshset_info
192  );
193  int meshset_data[2];
194  EntityHandle root = moab.get_root_set();
195  rval = moab.tag_get_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
196  if(meshset_data[0]==0) {
197  meshset_data[0] = 1;
198  rval = moab.tag_set_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
199 
200  }
201  bit_levels.push_back(BitRefLevel().set(meshset_data[0]-1));
202  }
203  BitRefLevel problem_bit_level = bit_levels.back();
204  ierr = m_field.seed_ref_level_3D(0,problem_bit_level); CHKERRQ(ierr);
205 
206  // const clock_t begin_time = clock();
207  ierr = m_field.build_fields(); CHKERRQ(ierr);
208  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
209 
210  Range preriodic_prisms;
211  if(choise_value == HOMOBCPERIODIC) {
212  //FIXME: Naming convention is not consistent in this section of code
213  //=======================================================================================================
214  //Add Periodic Prisims Between Triangles on -ve and +ve faces to implement periodic bounary conditions
215  //=======================================================================================================
216  //Populating the Multi-index container with -ve triangles
217  Range SurTrisNeg;
218  ierr = m_field.get_cubit_msId_entities_by_dimension(101,SIDESET,2,SurTrisNeg,true); CHKERRQ(ierr);
219  ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
220  Face_CenPos_Handle_multiIndex Face_CenPos_Handle_varNeg, Face_CenPos_Handle_varPos;
221  double TriCen[3], coords_Tri[9];
222 
223 
224  double roundfact=1e3;
225 
226  for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
227  const EntityHandle* conn_face; int num_nodes_Tri;
228 
229  //get nodes attached to the face
230  rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
231  //get nodal coordinates
232  rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
233 
234  //Find triangle centriod
235  TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
236  TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
237  TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
238 
239  //round values to roundfact disimal places
240  TriCen[0]=round(TriCen[0]*roundfact)/roundfact;
241  TriCen[1]=round(TriCen[1]*roundfact)/roundfact;
242  TriCen[2]=round(TriCen[2]*roundfact)/roundfact;
243 
244  // cout<<"\n\n\nTriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
245  //fill the multi-index container with centriod coordinates and triangle handles
246  Face_CenPos_Handle_varNeg.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
247  // for(int ii=0; ii<3; ii++) cout<<"TriCen "<<TriCen[ii]<<endl;
248  // cout<<endl<<endl;
249  }
250 
251  //Populating the Multi-index container with +ve triangles
252  Range SurTrisPos;
253  ierr = m_field.get_cubit_msId_entities_by_dimension(102,SIDESET,2,SurTrisPos,true); CHKERRQ(ierr);
254  ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 102 = %d\n",SurTrisPos.size()); CHKERRQ(ierr);
255  for(Range::iterator it = SurTrisPos.begin(); it!=SurTrisPos.end(); it++) {
256  const EntityHandle* conn_face; int num_nodes_Tri;
257 
258  //get nodes attached to the face
259  rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
260  //get nodal coordinates
261  rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
262 
263  //Find triangle centriod
264  TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
265  TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
266  TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
267 
268  //round values to roundfact disimal places
269  TriCen[0]=round(TriCen[0]*roundfact)/roundfact;
270  TriCen[1]=round(TriCen[1]*roundfact)/roundfact;
271  TriCen[2]=round(TriCen[2]*roundfact)/roundfact;
272  // cout<<"\n\n\nTriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
273 
274  //fill the multi-index container with centriod coordinates and triangle handles
275  Face_CenPos_Handle_varPos.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
276  }
277 
278  //Find minimum and maximum X, Y and Z coordinates of the RVE (using multi-index container)
279  double XcoordMin, YcoordMin, ZcoordMin, XcoordMax, YcoordMax, ZcoordMax;
280  typedef Face_CenPos_Handle_multiIndex::index<xcoord_tag>::type::iterator Tri_Xcoord_iterator;
281  typedef Face_CenPos_Handle_multiIndex::index<ycoord_tag>::type::iterator Tri_Ycoord_iterator;
282  typedef Face_CenPos_Handle_multiIndex::index<zcoord_tag>::type::iterator Tri_Zcoord_iterator;
283  Tri_Xcoord_iterator XcoordMin_it, XcoordMax_it;
284  Tri_Ycoord_iterator YcoordMin_it, YcoordMax_it;
285  Tri_Zcoord_iterator ZcoordMin_it, ZcoordMax_it;
286 
287  //XcoordMax_it-- because .end() will point iterator after the data range but .begin() will point the iteratore to the first value of range
288  XcoordMin_it=Face_CenPos_Handle_varNeg.get<xcoord_tag>().begin(); XcoordMin=XcoordMin_it->xcoord;
289  XcoordMax_it=Face_CenPos_Handle_varPos.get<xcoord_tag>().end(); XcoordMax_it--; XcoordMax=XcoordMax_it->xcoord;
290  YcoordMin_it=Face_CenPos_Handle_varNeg.get<ycoord_tag>().begin(); YcoordMin=YcoordMin_it->ycoord;
291  YcoordMax_it=Face_CenPos_Handle_varPos.get<ycoord_tag>().end(); YcoordMax_it--; YcoordMax=YcoordMax_it->ycoord;
292  ZcoordMin_it=Face_CenPos_Handle_varNeg.get<zcoord_tag>().begin(); ZcoordMin=ZcoordMin_it->zcoord;
293  ZcoordMax_it=Face_CenPos_Handle_varPos.get<zcoord_tag>().end(); ZcoordMax_it--; ZcoordMax=ZcoordMax_it->zcoord;
294 
295  /*double LxRVE, LyRVE, LzRVE;
296  LxRVE=XcoordMax-XcoordMin;
297  LyRVE=YcoordMax-YcoordMin;
298  LzRVE=ZcoordMax-ZcoordMin;*/
299 
300  //Creating Prisims between triangles on -ve and +ve faces
301  typedef Face_CenPos_Handle_multiIndex::index<Tri_Hand_tag>::type::iterator Tri_Hand_iterator;
302  Tri_Hand_iterator Tri_Neg;
303  typedef Face_CenPos_Handle_multiIndex::index<Composite_xyzcoord>::type::iterator xyzcoord_iterator;
304  xyzcoord_iterator Tri_Pos;
305  double XPos, YPos, ZPos;
306  //int count=0;
307 
308  //loop over -ve triangles to create prisims elemenet between +ve and -ve triangles
309  for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
310 
311  Tri_Neg=Face_CenPos_Handle_varNeg.get<Tri_Hand_tag>().find(*it);
312  // cout<<"xcoord= "<<Tri_iit->xcoord << " ycoord= "<< Tri_iit->ycoord << " ycoord= "<< Tri_iit->zcoord <<endl;
313 
314  //corresponding +ve triangle
315  if(Tri_Neg->xcoord==XcoordMin){XPos=XcoordMax; YPos=Tri_Neg->ycoord; ZPos=Tri_Neg->zcoord;};
316  if(Tri_Neg->ycoord==YcoordMin){XPos=YPos=Tri_Neg->xcoord; YPos=YcoordMax; ZPos=Tri_Neg->zcoord;};
317  if(Tri_Neg->zcoord==ZcoordMin){XPos=YPos=Tri_Neg->xcoord; YPos=Tri_Neg->ycoord; ZPos=ZcoordMax; };
318  Tri_Pos=Face_CenPos_Handle_varPos.get<Composite_xyzcoord>().find(boost::make_tuple(XPos, YPos, ZPos));
319  // cout<<"Tri_Neg->xcoord= "<<Tri_Neg->xcoord << " Tri_Neg->ycoord "<< Tri_Neg->ycoord << " Tri_Neg->zcoord= "<< Tri_Neg->zcoord <<endl;
320  // cout<<"Tri_Pos->xcoord= "<<Tri_Pos->xcoord << " Tri_Pos->ycoord "<< Tri_Pos->ycoord << " Tri_Pos->zcoord= "<< Tri_Pos->zcoord <<endl;
321 
322 
323  //+ve and -ve nodes and their coords (+ve and -ve tiangles nodes can have matching problems, which can produce twisted prism)
324  EntityHandle PrismNodes[6];
325  vector<EntityHandle> TriNodesNeg, TriNodesPos;
326  double CoordNodeNeg[9], CoordNodePos[9];
327  rval = moab.get_connectivity(&(Tri_Neg->Tri_Hand),1,TriNodesNeg,true); CHKERRQ_MOAB(rval);
328  rval = moab.get_connectivity(&(Tri_Pos->Tri_Hand),1,TriNodesPos,true); CHKERRQ_MOAB(rval);
329  rval = moab.get_coords(&TriNodesNeg[0],3,CoordNodeNeg); MOAB_THROW(rval);
330  rval = moab.get_coords(&TriNodesPos[0],3,CoordNodePos); MOAB_THROW(rval);
331  for(int ii=0; ii<3; ii++){
332  PrismNodes[ii]=TriNodesNeg[ii];
333  }
334  // for(int ii=0; ii<3; ii++){
335  // cout<<"xcoord= "<<CoordNodeNeg[3*ii] << " ycoord= "<< CoordNodeNeg[3*ii+1] << " zcoord= "<< CoordNodeNeg[3*ii+2] <<endl;
336  // }
337  // for(int ii=0; ii<3; ii++){
338  // cout<<"xcoord= "<<CoordNodePos[3*ii] << " ycoord= "<< CoordNodePos[3*ii+1] << " zcoord= "<< CoordNodePos[3*ii+2] <<endl;
339  // }
340 
341  //Match exact nodes to each other to avoide the problem of twisted prisms
342  double XNodeNeg, YNodeNeg, ZNodeNeg, XNodePos, YNodePos, ZNodePos;
343  for(int ii=0; ii<3; ii++){
344  if(Tri_Neg->xcoord==XcoordMin){XNodeNeg=XcoordMax; YNodeNeg=CoordNodeNeg[3*ii+1]; ZNodeNeg=CoordNodeNeg[3*ii+2];};
345  if(Tri_Neg->ycoord==YcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=YcoordMax; ZNodeNeg=CoordNodeNeg[3*ii+2];};
346  if(Tri_Neg->zcoord==ZcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=CoordNodeNeg[3*ii+1]; ZNodeNeg=ZcoordMax;};
347  for(int jj=0; jj<3; jj++){
348  //Round nodal coordinates to 3 dicimal places only for comparison //round values to 3 disimal places
349  XNodeNeg=round(XNodeNeg*roundfact)/roundfact;
350  YNodeNeg=round(YNodeNeg*roundfact)/roundfact;
351  ZNodeNeg=round(ZNodeNeg*roundfact)/roundfact;
352 
353 
354  XNodePos=CoordNodePos[3*jj]; YNodePos=CoordNodePos[3*jj+1]; ZNodePos=CoordNodePos[3*jj+2];
355  XNodePos=round(XNodePos*roundfact)/roundfact;
356  YNodePos=round(YNodePos*roundfact)/roundfact;
357  ZNodePos=round(ZNodePos*roundfact)/roundfact;
358 
359  if(XNodeNeg==XNodePos && YNodeNeg==YNodePos && ZNodeNeg==ZNodePos){
360  PrismNodes[3+ii]=TriNodesPos[jj];
361  break;
362  }
363  }
364  }
365  //prism nodes and their coordinates
366  double CoordNodesPrisms[18];
367  rval = moab.get_coords(PrismNodes,6,CoordNodesPrisms); MOAB_THROW(rval);
368  // for(int ii=0; ii<6; ii++){
369  // cout<<"xcoord= "<<CoordNodesPrisms[3*ii] << " ycoord= "<< CoordNodesPrisms[3*ii+1] << " zcoord= "<< CoordNodesPrisms[3*ii+2] <<endl;
370  // }
371  // cout<<endl<<endl;
372  //insertion of individula prism element and its addition to range preriodic_prisms
373  EntityHandle PeriodicPrism;
374  rval = moab.create_element(MBPRISM,PrismNodes,6,PeriodicPrism); CHKERRQ_MOAB(rval);
375  preriodic_prisms.insert(PeriodicPrism);
376 
377  }
378 
379  //insertion of individual prism element and its addition to range preriodic_prisms
380  ierr = m_field.seed_ref_level(preriodic_prisms,problem_bit_level); CHKERRQ(ierr);
381 
382  }
383 
384  EntityHandle out_meshset;
385  rval = moab.create_meshset(MESHSET_SET,out_meshset); CHKERRQ_MOAB(rval);
386  // ierr = m_field.get_problem_finite_elements_entities("POTENTIAL_PROBLEM","POTENTIAL_ELEM",out_meshset); CHKERRQ(ierr);
387  ierr = m_field.get_entities_by_ref_level(bit_levels.back(),BitRefLevel().set(),out_meshset); CHKERRQ(ierr);
388  Range LatestRefinedTets;
389  rval = moab.get_entities_by_type(out_meshset, MBTET,LatestRefinedTets,true); CHKERRQ_MOAB(rval);
390  Range LatestRefinedPrisms;
391  rval = moab.get_entities_by_type(out_meshset, MBPRISM,LatestRefinedPrisms,true); CHKERRQ_MOAB(rval);
392 
393  Range prims_on_problem_bit_level;
394  ierr = m_field.get_entities_by_type_and_ref_level(
395  problem_bit_level,BitRefLevel().set(),MBPRISM,prims_on_problem_bit_level
396  ); CHKERRQ(ierr);
397  //to create meshset from range
398  EntityHandle meshset_prims_on_problem_bit_level;
399  rval = moab.create_meshset(MESHSET_SET,meshset_prims_on_problem_bit_level); CHKERRQ_MOAB(rval);
400  rval = moab.add_entities(meshset_prims_on_problem_bit_level,prims_on_problem_bit_level); CHKERRQ_MOAB(rval);
401  ierr = m_field.seed_ref_level_MESHSET(meshset_prims_on_problem_bit_level,BitRefLevel().set()); CHKERRQ(ierr);
402 
403  //Fields
404  int field_rank=3;
405  ierr = m_field.add_field("DISPLACEMENT",H1,AINSWORTH_LEGENDRE_BASE,field_rank); CHKERRQ(ierr);
406  ierr = m_field.add_ents_to_field_by_type(0,MBTET,"DISPLACEMENT"); CHKERRQ(ierr);
407  ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,field_rank,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
408  ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
409 
410  //add entitities (by tets) to the field
411  if(choise_value == HOMOBCDISP) {
412  ierr = m_field.add_field("LAGRANGE_MUL_DISP",H1,AINSWORTH_LEGENDRE_BASE,field_rank); CHKERRQ(ierr);
413  for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
414  if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
415  Range tris;
416  rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
417  ierr = m_field.add_ents_to_field_by_type(tris,MBTRI,"LAGRANGE_MUL_DISP"); CHKERRQ(ierr);
418  }
419  }
420  }
421 
422  if(choise_value == HOMOBCPERIODIC) {
423  ierr = m_field.add_field("LAGRANGE_MUL_PERIODIC",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
424  //Control 3 rigid body translations in x, y and z axis
425  ierr = m_field.add_field("LAGRANGE_MUL_RIGID_TRANS",NOFIELD,NOBASE,3); CHKERRQ(ierr);
426  // Setting up dummy no-field vertex
427  EntityHandle no_field_vertex;
428  {
429  const double coords[] = {0,0,0};
430  rval = m_field.get_moab().create_vertex(coords,no_field_vertex); CHKERRQ_MOAB(rval);
431  Range range_no_field_vertex;
432  range_no_field_vertex.insert(no_field_vertex);
433  ierr = m_field.seed_ref_level(range_no_field_vertex,BitRefLevel().set()); CHKERRQ(ierr);
434  EntityHandle meshset;
435  meshset = m_field.get_field_meshset("LAGRANGE_MUL_RIGID_TRANS");
436  rval = m_field.get_moab().add_entities(meshset,range_no_field_vertex); CHKERRQ_MOAB(rval);
437  }
438  Range surface_negative;
439  ierr = m_field.get_cubit_msId_entities_by_dimension(101,SIDESET,2,surface_negative,true); CHKERRQ(ierr);
440  ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 101 = %d\n",surface_negative.size()); CHKERRQ(ierr);
441  ierr = m_field.add_ents_to_field_by_type(surface_negative,MBTRI,"LAGRANGE_MUL_PERIODIC"); CHKERRQ(ierr);
442  }
443 
444  if(choise_value == HOMOBCTRAC) {
445  ierr = m_field.add_field("LAGRANGE_MUL_TRAC",NOFIELD,NOBASE,6); CHKERRQ(ierr);
446  //Control 3 rigid body translations in x, y and z axis
447  ierr = m_field.add_field("LAGRANGE_MUL_RIGID_TRANS",NOFIELD,NOBASE,3); CHKERRQ(ierr);
448  //Controla 3 rigid body rotations about x, y and z axis
449  ierr = m_field.add_field("LAGRANGE_MUL_RIGID_ROT",NOFIELD,NOBASE,3); CHKERRQ(ierr);
450  EntityHandle no_field_vertex;
451  {
452  const double coords[] = {0,0,0};
453  rval = m_field.get_moab().create_vertex(coords,no_field_vertex); CHKERRQ_MOAB(rval);
454  Range range_no_field_vertex;
455  range_no_field_vertex.insert(no_field_vertex);
456  ierr = m_field.seed_ref_level(range_no_field_vertex,BitRefLevel().set()); CHKERRQ(ierr);
457  EntityHandle meshset;
458  meshset = m_field.get_field_meshset("LAGRANGE_MUL_TRAC");
459  rval = m_field.get_moab().add_entities(meshset,range_no_field_vertex); CHKERRQ_MOAB(rval);
460  meshset = m_field.get_field_meshset("LAGRANGE_MUL_RIGID_TRANS");
461  rval = m_field.get_moab().add_entities(meshset,range_no_field_vertex); CHKERRQ_MOAB(rval);
462  meshset = m_field.get_field_meshset("LAGRANGE_MUL_RIGID_ROT");
463  rval = m_field.get_moab().add_entities(meshset,range_no_field_vertex); CHKERRQ_MOAB(rval);
464  }
465  }
466 
467  //set app. order
468  //see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes (Mark Ainsworth & Joe Coyle)
469  ierr = m_field.set_field_order(0,MBTET,"DISPLACEMENT",order); CHKERRQ(ierr);
470  ierr = m_field.set_field_order(0,MBTRI,"DISPLACEMENT",order); CHKERRQ(ierr);
471  ierr = m_field.set_field_order(0,MBEDGE,"DISPLACEMENT",order); CHKERRQ(ierr);
472  ierr = m_field.set_field_order(0,MBVERTEX,"DISPLACEMENT",1); CHKERRQ(ierr);
473 
474  ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
475  ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
476  ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
477  ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
478 
479  PetscBool fo_boundary = PETSC_FALSE;
480  ierr = PetscOptionsGetBool(PETSC_NULL,"-my_fo_boundary",&fo_boundary,PETSC_NULL); CHKERRQ(ierr);
481  if(fo_boundary) {
482  for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
483  if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
484  Range tris;
485  rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
486  Range tris_edges;
487  rval = moab.get_adjacencies(tris,1,false,tris_edges,moab::Interface::UNION); CHKERRQ_MOAB(rval);
488  ierr = m_field.set_field_order(tris,"DISPLACEMENT",1); CHKERRQ(ierr);
489  ierr = m_field.set_field_order(tris_edges,"DISPLACEMENT",1); CHKERRQ(ierr);
490  }
491  }
492  }
493 
494  if(choise_value == HOMOBCDISP) {
495  ierr = m_field.set_field_order(0,MBTRI,"LAGRANGE_MUL_DISP",order); CHKERRQ(ierr);
496  ierr = m_field.set_field_order(0,MBEDGE,"LAGRANGE_MUL_DISP",order); CHKERRQ(ierr);
497  ierr = m_field.set_field_order(0,MBVERTEX,"LAGRANGE_MUL_DISP",1); CHKERRQ(ierr);
498  }
499 
500  if(choise_value == HOMOBCPERIODIC) {
501  ierr = m_field.set_field_order(0,MBTRI,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
502  ierr = m_field.set_field_order(0,MBEDGE,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
503  ierr = m_field.set_field_order(0,MBVERTEX,"LAGRANGE_MUL_PERIODIC",1); CHKERRQ(ierr);
504  }
505 
506  //build field
507  ierr = m_field.build_fields(); CHKERRQ(ierr);
508 
509  Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
510  ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
511 
512  boost::shared_ptr<Hooke<adouble> > hooke_adouble(new Hooke<adouble>());
513  boost::shared_ptr<Hooke<double> > hooke_double(new Hooke<double>());
514 
515  NonlinearElasticElement iso_elastic(m_field,1);
516  for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,BLOCKSET|MAT_ELASTICSET,it)) {
517  if(it->getName() != "MAT_ELASTIC_1") continue;
518  Mat_Elastic mydata;
519  ierr = it->getAttributeDataStructure(mydata); CHKERRQ(ierr);
520  int id = it->getMeshsetId();
521  EntityHandle meshset = it->getMeshset();
522  rval = m_field.get_moab().get_entities_by_type(
523  meshset,MBTET,iso_elastic.setOfBlocks[id].tEts,true
524  ); CHKERRQ_MOAB(rval);
525  iso_elastic.setOfBlocks[id].iD = id;
526  iso_elastic.setOfBlocks[id].E = mydata.data.Young;
527  iso_elastic.setOfBlocks[id].PoissonRatio = mydata.data.Poisson;
528  iso_elastic.setOfBlocks[id].materialDoublePtr = hooke_double;
529  iso_elastic.setOfBlocks[id].materialAdoublePtr = hooke_adouble;
530  ierr = m_field.seed_finite_elements(iso_elastic.setOfBlocks[id].tEts); CHKERRQ(ierr);
531  }
532  ierr = iso_elastic.addElement("ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
533  ierr = iso_elastic.setOperators("DISPLACEMENT","MESH_NODE_POSITIONS",false,true); CHKERRQ(ierr);
534  if(m_field.check_field("POTENTIAL_FIELD")) {
535  ierr = m_field.modify_finite_element_add_field_data("ELASTIC","POTENTIAL_FIELD"); CHKERRQ(ierr);
536  }
537 
538  NonlinearElasticElement trans_elastic(m_field,2);
539  trans_elastic.commonData.spatialPositions = "DISPLACEMENT";
540  trans_elastic.commonData.meshPositions = "MESH_NODE_POSITIONS";
541  std::map<int,boost::shared_ptr<SmallStrainTranverslyIsotropicADouble> > tranversly_isotropic_adouble_ptr_map;
542  std::map<int,boost::shared_ptr<SmallStrainTranverslyIsotropicDouble> > tranversly_isotropic_double_ptr_map;
543  bool trans_iso_blocks = false;
544  for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
545  //Get block name
546  string name = it->getName();
547  if (name.compare(0,20,"MAT_ELASTIC_TRANSISO") == 0) {
548  trans_iso_blocks = true;
549  int id = it->getMeshsetId();
550  Mat_Elastic_TransIso mydata;
551  ierr = it->getAttributeDataStructure(mydata); CHKERRQ(ierr);
552  tranversly_isotropic_adouble_ptr_map[id] = boost::make_shared<SmallStrainTranverslyIsotropicADouble>();
553  tranversly_isotropic_double_ptr_map[id] = boost::make_shared<SmallStrainTranverslyIsotropicDouble>();
554  //nu_p, nu_pz, E_p, E_z, G_zp
555  tranversly_isotropic_adouble_ptr_map.at(id)->E_p = mydata.data.Youngp;
556  tranversly_isotropic_double_ptr_map.at(id)->E_p = mydata.data.Youngp;
557  tranversly_isotropic_adouble_ptr_map.at(id)->E_z = mydata.data.Youngz;
558  tranversly_isotropic_double_ptr_map.at(id)->E_z = mydata.data.Youngz;
559  tranversly_isotropic_adouble_ptr_map.at(id)->nu_p = mydata.data.Poissonp;
560  tranversly_isotropic_double_ptr_map.at(id)->nu_p = mydata.data.Poissonp;
561  tranversly_isotropic_adouble_ptr_map.at(id)->nu_pz = mydata.data.Poissonpz;
562  tranversly_isotropic_double_ptr_map.at(id)->nu_pz = mydata.data.Poissonpz;
563  double shear_zp;
564  if(mydata.data.Shearzp!=0) {
565  shear_zp = mydata.data.Shearzp;
566  } else {
567  shear_zp = mydata.data.Youngz/(2*(1+mydata.data.Poissonpz));
568  }
569  tranversly_isotropic_adouble_ptr_map.at(it->getMeshsetId())->G_zp = shear_zp;
570  tranversly_isotropic_double_ptr_map.at(it->getMeshsetId())->G_zp = shear_zp;
571  //get tets from block where material is defined
572  EntityHandle meshset = it->getMeshset();
573  rval = m_field.get_moab().get_entities_by_type(
574  meshset,MBTET,trans_elastic.setOfBlocks[id].tEts,true
575  ); CHKERRQ_MOAB(rval);
576  //adding material to nonlinear class
577  trans_elastic.setOfBlocks[id].iD = id;
578  //note that material parameters are defined internally in material model
579  trans_elastic.setOfBlocks[id].E = 0; // this is not working for this material
580  trans_elastic.setOfBlocks[id].PoissonRatio = 0; // this is not working for this material
581  trans_elastic.setOfBlocks[id].materialDoublePtr = tranversly_isotropic_double_ptr_map.at(id);
582  trans_elastic.setOfBlocks[id].materialAdoublePtr = tranversly_isotropic_adouble_ptr_map.at(id);
583  ierr = m_field.seed_finite_elements(trans_elastic.setOfBlocks[id].tEts); CHKERRQ(ierr);
584  }
585  }
586  if(trans_iso_blocks) {
587  ierr = m_field.add_finite_element("TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
588  ierr = m_field.add_finite_element("TRAN_ISOTROPIC_ELASTIC",MF_ZERO); CHKERRQ(ierr);
589  ierr = m_field.modify_finite_element_add_field_row("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
590  ierr = m_field.modify_finite_element_add_field_col("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
591  ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
592  ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","POTENTIAL_FIELD"); CHKERRQ(ierr);
593  if(m_field.check_field("MESH_NODE_POSITIONS")) {
594  ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
595  }
596  for(
597  map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
598  sit!=trans_elastic.setOfBlocks.end();sit++
599  ) {
600  ierr = m_field.add_ents_to_finite_element_by_type(sit->second.tEts,MBTET,"TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
601  }
602  }
603  if(trans_iso_blocks) {
604  //Rhs
605  trans_elastic.feRhs.getOpPtrVector().push_back(
606  new NonlinearElasticElement::OpGetCommonDataAtGaussPts("DISPLACEMENT",trans_elastic.commonData)
607  );
608  trans_elastic.feRhs.getOpPtrVector().push_back(
609  new NonlinearElasticElement::OpGetCommonDataAtGaussPts("POTENTIAL_FIELD",trans_elastic.commonData)
610  );
611  if(m_field.check_field("MESH_NODE_POSITIONS")) {
612  trans_elastic.feRhs.getOpPtrVector().push_back(
613  new NonlinearElasticElement::OpGetCommonDataAtGaussPts("MESH_NODE_POSITIONS",trans_elastic.commonData)
614  );
615  }
616  map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
617  for(;sit!=trans_elastic.setOfBlocks.end();sit++) {
618  trans_elastic.feRhs.getOpPtrVector().push_back(
619  new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
620  "DISPLACEMENT",sit->second,trans_elastic.commonData,2,false,false,true
621  )
622  );
623  trans_elastic.feRhs.getOpPtrVector().push_back(
624  new NonlinearElasticElement::OpRhsPiolaKirchhoff(
625  "DISPLACEMENT",sit->second,trans_elastic.commonData
626  )
627  );
628  }
629 
630  //Lhs
631  trans_elastic.feLhs.getOpPtrVector().push_back(
632  new NonlinearElasticElement::OpGetCommonDataAtGaussPts("DISPLACEMENT",trans_elastic.commonData)
633  );
634  trans_elastic.feLhs.getOpPtrVector().push_back(
635  new NonlinearElasticElement::OpGetCommonDataAtGaussPts("POTENTIAL_FIELD",trans_elastic.commonData)
636  );
637  if(m_field.check_field("MESH_NODE_POSITIONS")) {
638  trans_elastic.feLhs.getOpPtrVector().push_back(
639  new NonlinearElasticElement::OpGetCommonDataAtGaussPts("MESH_NODE_POSITIONS",trans_elastic.commonData)
640  );
641  }
642  sit = trans_elastic.setOfBlocks.begin();
643  for(;sit!=trans_elastic.setOfBlocks.end();sit++) {
644  trans_elastic.feLhs.getOpPtrVector().push_back(
645  new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
646  "DISPLACEMENT",sit->second,trans_elastic.commonData,2,true,false,true
647  )
648  );
649  trans_elastic.feLhs.getOpPtrVector().push_back(
650  new NonlinearElasticElement::OpLhsPiolaKirchhoff_dx(
651  "DISPLACEMENT","DISPLACEMENT",sit->second,trans_elastic.commonData
652  )
653  );
654  }
655  }
656 
657  BCs_RVELagrange_Disp lagrangian_element_disp(m_field);
658  if(choise_value == HOMOBCDISP) {
659  lagrangian_element_disp.addLagrangiangElement(
660  "LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS"
661  );
662  }
663 
664  BCs_RVELagrange_Trac lagrangian_element_trac(m_field);
665  BCs_RVELagrange_Trac_Rigid_Trans lagrangian_element_rigid_body_trans(m_field);
666  BCs_RVELagrange_Trac_Rigid_Rot lagrangian_element_rigid_body_rot(m_field);
667  if(choise_value == HOMOBCTRAC) {
668  lagrangian_element_trac.addLagrangiangElement(
669  "LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_TRAC","MESH_NODE_POSITIONS"
670  );
671  lagrangian_element_trac.addLagrangiangElement(
672  "LAGRANGE_ELEM_TRANS","DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS","MESH_NODE_POSITIONS"
673  );
674  lagrangian_element_trac.addLagrangiangElement(
675  "LAGRANGE_ELEM_ROT","DISPLACEMENT","LAGRANGE_MUL_RIGID_ROT","MESH_NODE_POSITIONS"
676  );
677  }
678 
679  BCs_RVELagrange_Periodic lagrangian_element_periodic(m_field);
680  if(choise_value == HOMOBCPERIODIC) {
681  lagrangian_element_periodic.addLagrangiangElement(
682  "LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",preriodic_prisms
683  );
684  lagrangian_element_trac.addLagrangiangElement(
685  "LAGRANGE_ELEM_TRANS","DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS","MESH_NODE_POSITIONS"
686  );
687  }
688 
689  struct MinMaxNodes {
690  enum MINAMX { C0,MAXLAST };
691  EntityHandle entMinMax[MAXLAST];
692  ublas::vector<int> rowIndices;
693  VectorDouble rHs[6];
694  MinMaxNodes() {
695  rowIndices.resize(3*MAXLAST);
696  for(int ii = 0;ii<6;ii++) {
697  rHs[ii].resize(3*MAXLAST);
698  }
699  }
700  };
701  MinMaxNodes minMaxNodes;
702 
703  if(choise_value == NITSCHE) { // Condensed traction/periodc BC
704  ierr = m_field.add_finite_element("SURFACE_ELEMENTS"); CHKERRQ(ierr);
705  ierr = m_field.modify_finite_element_add_field_row("SURFACE_ELEMENTS","DISPLACEMENT"); CHKERRQ(ierr);
706  ierr = m_field.modify_finite_element_add_field_col("SURFACE_ELEMENTS","DISPLACEMENT"); CHKERRQ(ierr);
707  ierr = m_field.modify_finite_element_add_field_data("SURFACE_ELEMENTS","DISPLACEMENT"); CHKERRQ(ierr);
708  ierr = m_field.modify_finite_element_add_field_data("SURFACE_ELEMENTS","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
709  EntityHandle condensed_traction_element_meshset;
710  rval = moab.create_meshset(MESHSET_TRACK_OWNER,condensed_traction_element_meshset); CHKERRQ(ierr);
711  Range nodes;
712  for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
713  if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
714  Range tris;
715  rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
716  ierr = m_field.add_ents_to_finite_element_by_type(tris,MBTRI,"SURFACE_ELEMENTS"); CHKERRQ(ierr);
717  Range tris_nodes;
718  rval = moab.get_connectivity(tris,nodes,true); CHKERRQ_MOAB(rval);
719  nodes.merge(tris_nodes);
720  }
721  }
722 
723  {
724  VectorDouble x,y,z;
725  x.resize(nodes.size(),false);
726  y.resize(nodes.size(),false);
727  z.resize(nodes.size(),false);
728  rval = moab.get_coords(nodes,&x[0],&y[0],&z[0]); CHKERRQ_MOAB(rval);
729  int bc_nb = 0;
730  for(int sx = -1; sx<=+1; sx+=2) {
731  for(int sy = -1; sy<=+1; sy+=2) {
732  for(int sz = -1; sz<=+1; sz+=2) {
733  if(bc_nb == MinMaxNodes::MAXLAST) break;
734  VectorDouble dist_up_right;
735  dist_up_right.resize(x.size(),false);
736  dist_up_right.clear();
737  for(unsigned int nn = 0;nn<x.size();nn++) {
738  if(
739  ((sx*x[nn])>0)&&
740  ((sy*y[nn])>0)&&
741  ((sz*z[nn])>0)
742  ) {
743  dist_up_right[nn] = sx*x[nn]+sy*y[nn]+sz*z[nn];
744  }
745  }
746  VectorDouble::iterator dist_it;
747  dist_it = max_element(dist_up_right.begin(),dist_up_right.end());
748  minMaxNodes.entMinMax[bc_nb++] = nodes[std::distance(dist_up_right.begin(),dist_it)];
749  }
750  }
751  }
752  }
753 
754  }
755 
756  //build finite elements
757  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
758  //build adjacencies
759  ierr = m_field.build_adjacencies(problem_bit_level); CHKERRQ(ierr);
760 
761  //define problems
762  ierr = m_field.add_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
763  //set finite elements for problem
764  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","ELASTIC"); CHKERRQ(ierr);
765  if(trans_iso_blocks) {
766  ierr = m_field.modify_problem_add_finite_element(
767  "ELASTIC_MECHANICS","TRAN_ISOTROPIC_ELASTIC"
768  ); CHKERRQ(ierr);
769  }
770  if(choise_value == HOMOBCDISP) {
771  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM"); CHKERRQ(ierr);
772  }
773  if(choise_value == HOMOBCTRAC) {
774  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM"); CHKERRQ(ierr);
775  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
776  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM_ROT"); CHKERRQ(ierr);
777  }
778  if(choise_value == HOMOBCPERIODIC) {
779  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM"); CHKERRQ(ierr);
780  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
781  }
782  if(choise_value == NITSCHE) {
783  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","SURFACE_ELEMENTS"); CHKERRQ(ierr);
784  }
785 
786  //set refinement level for problem
787  ierr = m_field.modify_problem_ref_level_add_bit("ELASTIC_MECHANICS",problem_bit_level); CHKERRQ(ierr);
788  //build problem
789  ierr = m_field.build_problems(); CHKERRQ(ierr);
790 
791  /****/
792  //mesh partitioning
793 
794  //partition
795  ierr = m_field.partition_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
796  ierr = m_field.partition_finite_elements(
797  "ELASTIC_MECHANICS",false,0,m_field.get_comm_size() // build elements on all procs
798  ); CHKERRQ(ierr);
799  //what are ghost nodes, see Petsc Manual
800  ierr = m_field.partition_ghost_dofs("ELASTIC_MECHANICS"); CHKERRQ(ierr);
801 
802  //create matrices
803  Vec D;
804  vector<Vec> F(6);
805  ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",ROW,&F[0]); CHKERRQ(ierr);
806  for(int ii = 1;ii<6;ii++) {
807  ierr = VecDuplicate(F[0],&F[ii]); CHKERRQ(ierr);
808  }
809  ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",COL,&D); CHKERRQ(ierr);
810 
811  Mat Aij;
812  ierr = m_field.MatCreateMPIAIJWithArrays("ELASTIC_MECHANICS",&Aij); CHKERRQ(ierr);
813  ierr = MatSetOption(Aij,MAT_STRUCTURALLY_SYMMETRIC,PETSC_TRUE); CHKERRQ(ierr);
814  ierr = MatSetOption(Aij,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE); CHKERRQ(ierr);
815  ierr = MatSetOption(Aij,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE); CHKERRQ(ierr);
816  ierr = MatSetOption(Aij,MAT_USE_INODES,PETSC_TRUE); CHKERRQ(ierr);
817  ierr = MatSetOption(Aij,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE); CHKERRQ(ierr);
818  ierr = MatSetOption(Aij,MAT_KEEP_NONZERO_PATTERN,PETSC_FALSE); CHKERRQ(ierr);
819 
820  /*{
821  ierr = MatView(Aij,PETSC_VIEWER_DRAW_SELF); CHKERRQ(ierr);
822  std::string wait;
823  std::cin >> wait;
824  }*/
825 
826  ierr = VecZeroEntries(D); CHKERRQ(ierr);
827  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
828  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
829  ierr = m_field.set_global_ghost_vector(
830  "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
831  ); CHKERRQ(ierr);
832  for(int ii = 0;ii<6;ii++) {
833  ierr = VecZeroEntries(F[ii]); CHKERRQ(ierr);
834  ierr = VecGhostUpdateBegin(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
835  ierr = VecGhostUpdateEnd(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
836  }
837  ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
838 
839  // NitscheMethod::BlockData nitsche_block_data;
840  // NitscheMethod::CommonData nitsche_common_data;
841  // PeriodicNitscheConstrains::CommonData periodic_common_data;
842  // PeriodicNitscheConstrains::MyNitscheVolume nitsche_element_iso(
843  // m_field,nitsche_block_data,nitsche_common_data,periodic_common_data
844  // );
845  // PeriodicNitscheConstrains::MyNitscheVolume nitsche_element_trans(
846  // m_field,nitsche_block_data,nitsche_common_data,periodic_common_data
847  // );
848 
849  // if(choise_value == NITSCHE) {
850 
851  // nitsche_block_data.faceElemName = "SURFACE_ELEMENTS";
852  // for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
853  // if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
854  // Range tris;
855  // rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
856  // nitsche_block_data.fAces.merge(tris);
857  // }
858  // }
859 
860  // for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
861  // if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
862  // Range tris;
863  // rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
864  // periodic_common_data.skinFaces.merge(tris);
865  // }
866  // }
867 
868  // nitsche_block_data.gamma = 1e-7;
869  // nitsche_block_data.phi = -1;
870  // periodic_common_data.eps = 0;
871  // ierr = PetscOptionsGetReal(
872  // PETSC_NULL,"-my_gamma",&nitsche_block_data.gamma,PETSC_NULL
873  // ); CHKERRQ(ierr);
874  // ierr = PetscOptionsGetReal(
875  // PETSC_NULL,"-my_phi",&nitsche_block_data.phi,PETSC_NULL
876  // ); CHKERRQ(ierr);
877  // ierr = PetscOptionsGetReal(
878  // PETSC_NULL,"-my_eps",&periodic_common_data.eps,PETSC_NULL
879  // ); CHKERRQ(ierr);
880  // ierr = PetscPrintf(
881  // PETSC_COMM_WORLD,
882  // "Nitsche method gamma = %4.2e phi = %2.1f eps = %4.2e\n",
883  // nitsche_block_data.gamma,nitsche_block_data.phi,periodic_common_data.eps
884  // ); CHKERRQ(ierr);
885 
886  // for(
887  // map<int,NonlinearElasticElement::BlockData>::iterator mit = iso_elastic.setOfBlocks.begin();
888  // mit!=iso_elastic.setOfBlocks.end();
889  // mit++
890  // ) {
891  // NonlinearElasticElement::CommonData &elastic_common_data = iso_elastic.commonData;
892  // NonlinearElasticElement::BlockData &elastic_block_data = mit->second;
893  // nitsche_element_iso.snes_B = Aij;
894 
895  // nitsche_element_iso.getOpPtrVector().push_back(
896  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
897  // "DISPLACEMENT",elastic_common_data
898  // )
899  // );
900  // nitsche_element_iso.getOpPtrVector().push_back(
901  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
902  // "MESH_NODE_POSITIONS",elastic_common_data
903  // )
904  // );
905  // nitsche_element_iso.getOpPtrVector().push_back(
906  // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
907  // "DISPLACEMENT",elastic_block_data,elastic_common_data,1,true,false,true
908  // )
909  // );
910  // nitsche_element_iso.getOpPtrVector().push_back(
911  // new PeriodicNitscheConstrains::OpLhsPeriodicNormal(
912  // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
913  // elastic_block_data,elastic_common_data,
914  // periodic_common_data
915  // )
916  // );
917  // nitsche_element_iso.getOpPtrVector().push_back(
918  // new PeriodicNitscheConstrains::OpRhsPeriodicNormal(
919  // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
920  // elastic_block_data,elastic_common_data,
921  // periodic_common_data,
922  // F
923  // )
924  // );
925  // nitsche_element_iso.getOpPtrVector().push_back(
926  // new NitscheMethod::OpLhsNormal(
927  // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
928  // elastic_block_data,elastic_common_data,true
929  // )
930  // );
931 
932  // // this is to get data on opposite element
933  // nitsche_element_iso.periodicVolume.getOpPtrVector().clear();
934  // nitsche_element_iso.periodicVolume.getOpPtrVector().push_back(
935  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
936  // "DISPLACEMENT",elastic_common_data
937  // )
938  // );
939  // nitsche_element_iso.periodicVolume.getOpPtrVector().push_back(
940  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
941  // "MESH_NODE_POSITIONS",elastic_common_data
942  // )
943  // );
944  // nitsche_element_iso.periodicVolume.getOpPtrVector().push_back(
945  // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
946  // "DISPLACEMENT",elastic_block_data,elastic_common_data,1,true,false,true
947  // )
948  // );
949  // nitsche_element_iso.periodicVolume.getOpPtrVector().push_back(
950  // new PeriodicNitscheConstrains::OpGetVolumeData(
951  // elastic_common_data,
952  // periodic_common_data
953  // )
954  // );
955  // periodic_common_data.volumeElemName = "ELASTIC";
956 
957  // nitsche_element_iso.addToRule = 1;
958  // ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC",nitsche_element_iso); CHKERRQ(ierr);
959  // }
960 
961  // for(
962  // map<int,NonlinearElasticElement::BlockData>::iterator mit = trans_elastic.setOfBlocks.begin();
963  // mit!=trans_elastic.setOfBlocks.end();
964  // mit++
965  // ) {
966  // NonlinearElasticElement::CommonData &elastic_common_data = trans_elastic.commonData;
967  // NonlinearElasticElement::BlockData &elastic_block_data = mit->second;
968  // nitsche_element_trans.snes_B = Aij;
969 
970  // nitsche_element_trans.getOpPtrVector().push_back(
971  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
972  // "DISPLACEMENT",elastic_common_data
973  // )
974  // );
975  // nitsche_element_trans.getOpPtrVector().push_back(
976  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
977  // "MESH_NODE_POSITIONS",elastic_common_data
978  // )
979  // );
980  // nitsche_element_trans.getOpPtrVector().push_back(
981  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
982  // "POTENTIAL_FIELD",elastic_common_data
983  // )
984  // );
985  // nitsche_element_trans.getOpPtrVector().push_back(
986  // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
987  // "DISPLACEMENT",elastic_block_data,elastic_common_data,2,true,false,true
988  // )
989  // );
990  // nitsche_element_trans.getOpPtrVector().push_back(
991  // new PeriodicNitscheConstrains::OpLhsPeriodicNormal(
992  // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
993  // elastic_block_data,elastic_common_data,
994  // periodic_common_data
995  // )
996  // );
997  // nitsche_element_trans.getOpPtrVector().push_back(
998  // new PeriodicNitscheConstrains::OpRhsPeriodicNormal(
999  // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
1000  // elastic_block_data,elastic_common_data,
1001  // periodic_common_data,
1002  // F
1003  // )
1004  // );
1005  // nitsche_element_trans.getOpPtrVector().push_back(
1006  // new NitscheMethod::OpLhsNormal(
1007  // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
1008  // elastic_block_data,elastic_common_data,true
1009  // )
1010  // );
1011 
1012  // // this is to get data on opposite element
1013  // nitsche_element_trans.periodicVolume.getOpPtrVector().clear();
1014  // nitsche_element_trans.periodicVolume.getOpPtrVector().push_back(
1015  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1016  // "DISPLACEMENT",elastic_common_data
1017  // )
1018  // );
1019  // nitsche_element_trans.periodicVolume.getOpPtrVector().push_back(
1020  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1021  // "POTENTIAL_FIELD",elastic_common_data
1022  // )
1023  // );
1024  // nitsche_element_trans.periodicVolume.getOpPtrVector().push_back(
1025  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1026  // "MESH_NODE_POSITIONS",elastic_common_data
1027  // )
1028  // );
1029  // nitsche_element_trans.periodicVolume.getOpPtrVector().push_back(
1030  // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
1031  // "DISPLACEMENT",elastic_block_data,elastic_common_data,2,true,false,true
1032  // )
1033  // );
1034  // nitsche_element_trans.periodicVolume.getOpPtrVector().push_back(
1035  // new PeriodicNitscheConstrains::OpGetVolumeData(
1036  // elastic_common_data,
1037  // periodic_common_data
1038  // )
1039  // );
1040  // periodic_common_data.volumeElemName = "TRAN_ISOTROPIC_ELASTIC";
1041 
1042  // nitsche_element_trans.addToRule = 1;
1043  // if(m_field.check_finite_element("TRAN_ISOTROPIC_ELASTIC"))
1044  // ierr = m_field.loop_finite_elements(
1045  // "ELASTIC_MECHANICS","TRAN_ISOTROPIC_ELASTIC",nitsche_element_trans
1046  // ); CHKERRQ(ierr);
1047  // }
1048 
1049 
1050  // }
1051 
1052  Vec volume_vec;
1053  int volume_vec_ghost[] = { 0 };
1054  ierr = VecCreateGhost(
1055  PETSC_COMM_WORLD,(!m_field.get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
1056  ); CHKERRQ(ierr);
1057  ierr = VecZeroEntries(volume_vec); CHKERRQ(ierr);
1058 
1059  iso_elastic.getLoopFeLhs().getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
1060  trans_elastic.getLoopFeLhs().getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
1061 
1062  //iso_elastic element matrix
1063  iso_elastic.getLoopFeLhs().snes_x = D;
1064  iso_elastic.getLoopFeLhs().snes_B = Aij;
1065  trans_elastic.getLoopFeLhs().snes_x = D;
1066  trans_elastic.getLoopFeLhs().snes_B = Aij;
1067  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC",iso_elastic.getLoopFeLhs()); CHKERRQ(ierr);
1068  if(m_field.check_finite_element("TRAN_ISOTROPIC_ELASTIC"))
1069  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","TRAN_ISOTROPIC_ELASTIC",trans_elastic.getLoopFeLhs()); CHKERRQ(ierr);
1070 
1071  if(choise_value == HOMOBCDISP) {
1072  lagrangian_element_disp.setRVEBCsOperators("DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS",Aij,F);
1073  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_disp.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
1074  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_disp.getLoopFeRVEBCRhs()); CHKERRQ(ierr);
1075  }
1076  if(choise_value == HOMOBCTRAC) {
1077  lagrangian_element_trac.setRVEBCsOperators("DISPLACEMENT","LAGRANGE_MUL_TRAC","MESH_NODE_POSITIONS",Aij,F);
1078  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_trac.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
1079  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_trac.getLoopFeRVEBCRhs()); CHKERRQ(ierr);
1080  lagrangian_element_rigid_body_trans.setRVEBCsRigidBodyTranOperators(
1081  "DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS",Aij,lagrangian_element_trac.setOfRVEBC
1082  );
1083  ierr = m_field.loop_finite_elements(
1084  "ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS",lagrangian_element_rigid_body_trans.getLoopFeRVEBCLhs()
1085  ); CHKERRQ(ierr);
1086  lagrangian_element_rigid_body_rot.setRVEBCsRigidBodyRotOperators(
1087  "DISPLACEMENT","LAGRANGE_MUL_RIGID_ROT",Aij,lagrangian_element_trac.setOfRVEBC
1088  );
1089  ierr = m_field.loop_finite_elements(
1090  "ELASTIC_MECHANICS","LAGRANGE_ELEM_ROT",lagrangian_element_rigid_body_rot.getLoopFeRVEBCLhs()
1091  ); CHKERRQ(ierr);
1092  }
1093  if(choise_value == HOMOBCPERIODIC) {
1094  lagrangian_element_periodic.setRVEBCsOperators("DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",Aij,F);
1095  ierr = m_field.loop_finite_elements(
1096  "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCLhs()
1097  ); CHKERRQ(ierr);
1098  ierr = m_field.loop_finite_elements(
1099  "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCRhs()
1100  ); CHKERRQ(ierr);
1101  lagrangian_element_rigid_body_trans.setRVEBCsRigidBodyTranOperators(
1102  "DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS",Aij,lagrangian_element_periodic.setOfRVEBC
1103  );
1104  ierr = m_field.loop_finite_elements(
1105  "ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS",lagrangian_element_rigid_body_trans.getLoopFeRVEBCLhs()
1106  ); CHKERRQ(ierr);
1107  }
1108 
1109  for(int ii = 0;ii<6;ii++) {
1110  ierr = VecGhostUpdateBegin(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
1111  ierr = VecGhostUpdateEnd(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
1112  ierr = VecAssemblyBegin(F[ii]); CHKERRQ(ierr);
1113  ierr = VecAssemblyEnd(F[ii]); CHKERRQ(ierr);
1114  }
1115  ierr = MatAssemblyBegin(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
1116  ierr = MatAssemblyEnd(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
1117 
1118  {
1119  //ierr = MatView(Aij,PETSC_VIEWER_DRAW_SELF); CHKERRQ(ierr);
1120  //std::string wait;
1121  //std::cin >> wait;
1122  }
1123 
1124  // if(choise_value == NITSCHE) {
1125  // if(periodic_common_data.eps==0) {
1126  // const Problem *problem_ptr;
1127  // ierr = m_field.get_problem("ELASTIC_MECHANICS",&problem_ptr); CHKERRQ(ierr);
1128  // for(int nn = 0;nn!=MinMaxNodes::MAXLAST;nn++) {
1129  // for(
1130  // _IT_NUMEREDDOF_ROW_BY_ENT_FOR_LOOP_(
1131  // problem_ptr,minMaxNodes.entMinMax[nn],dof
1132  // )
1133  // ) {
1134  // minMaxNodes.rowIndices[3*nn+dof->get()->getDofCoeffIdx()]
1135  // = dof->get()->getPetscGlobalDofIdx();
1136  // }
1137  // }
1138  // VectorDouble coords;
1139  // int nb_bcs = 3*MinMaxNodes::MAXLAST;
1140  // coords.resize(nb_bcs,false);
1141  // rval = moab.get_coords(&minMaxNodes.entMinMax[0],MinMaxNodes::MAXLAST,&coords[0]);
1142  // VectorDouble strain;
1143  // strain.resize(6,false);
1144  // MatrixDouble mat_d;
1145  // for(int ii = 0;ii<6;ii++) {
1146  // minMaxNodes.rHs[ii].clear();
1147  // strain.clear();
1148  // strain[ii] = 1;
1149  // for(int nn = 0;nn<MinMaxNodes::MAXLAST;nn++) {
1150  // PeriodicNitscheConstrains::OpRhsPeriodicNormal::calcualteDMatrix(
1151  // VectorAdaptor(3,ublas::shallow_array_adaptor<double>(3,&coords[3*nn])),
1152  // mat_d
1153  // );
1154  // VectorAdaptor rhs(3,ublas::shallow_array_adaptor<double>(3,&minMaxNodes.rHs[ii][3*nn]));
1155  // noalias(rhs) = prod(mat_d,strain);
1156  // }
1157  // }
1158  // for(int ii = 0;ii<6;ii++) {
1159  // ierr = VecSetValues(
1160  // F[ii],nb_bcs,&minMaxNodes.rowIndices[0],&minMaxNodes.rHs[ii][0],INSERT_VALUES
1161  // ); CHKERRQ(ierr);
1162  // ierr = VecAssemblyBegin(F[ii]); CHKERRQ(ierr);
1163  // ierr = VecAssemblyEnd(F[ii]); CHKERRQ(ierr);
1164  // }
1165  // ierr = MatZeroRows(
1166  // Aij,nb_bcs,&minMaxNodes.rowIndices[0],1,PETSC_NULL,PETSC_NULL
1167  // ); CHKERRQ(ierr);
1168  // }
1169  // }
1170 
1171  ierr = VecAssemblyBegin(volume_vec); CHKERRQ(ierr);
1172  ierr = VecAssemblyEnd(volume_vec); CHKERRQ(ierr);
1173  double rve_volume;
1174  ierr = VecSum(volume_vec,&rve_volume); CHKERRQ(ierr);
1175  ierr = PetscPrintf(PETSC_COMM_WORLD,"RVE Volume %3.2g\n",rve_volume); CHKERRQ(ierr);
1176 
1177  ierr = VecDestroy(&volume_vec);
1178 
1179  // Solver
1180  KSP solver;
1181  ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
1182  ierr = KSPSetOperators(solver,Aij,Aij); CHKERRQ(ierr);
1183  ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
1184  ierr = KSPSetUp(solver); CHKERRQ(ierr);
1185 
1186  MatrixDouble Dmat;
1187  Dmat.resize(6,6);
1188  Dmat.clear();
1189 
1190  //create a vector for 6 components of homogenized stress
1191  Vec stress_homo;
1192  int stress_homo_ghost[] = { 0,1,2,3,4,5,6 };
1193  NonlinearElasticElement::MyVolumeFE ave_stress_iso(m_field);
1194  NonlinearElasticElement::MyVolumeFE ave_stress_trans(m_field);
1195  PetscBool stress_by_boundary_integral = PETSC_FALSE;
1196  ierr = VecCreateGhost(
1197  PETSC_COMM_WORLD,(!m_field.get_comm_rank())?6:0,6,6,stress_homo_ghost,&stress_homo
1198  ); CHKERRQ(ierr);
1199  switch(choise_value) {
1200  case HOMOBCDISP:
1201  lagrangian_element_disp.setRVEBCsHomoStressOperators(
1202  "DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS",stress_homo
1203  );
1204  break;
1205  case HOMOBCTRAC:
1206  lagrangian_element_trac.setRVEBCsHomoStressOperators(
1207  "DISPLACEMENT","LAGRANGE_MUL_TRAC","MESH_NODE_POSITIONS",stress_homo
1208  );
1209  break;
1210  case HOMOBCPERIODIC:
1211  lagrangian_element_periodic.setRVEBCsHomoStressOperators(
1212  "DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",stress_homo
1213  );
1214  break;
1215  case NITSCHE:
1216  // if(stress_by_boundary_integral) {
1217  // nitsche_element_iso.getOpPtrVector().clear();
1218  // nitsche_element_trans.getOpPtrVector().clear();
1219  // nitsche_element_iso.periodicVolume.getOpPtrVector().clear();
1220  // nitsche_element_trans.periodicVolume.getOpPtrVector().clear();
1221  // for(
1222  // map<int,NonlinearElasticElement::BlockData>::iterator mit = iso_elastic.setOfBlocks.begin();
1223  // mit!=iso_elastic.setOfBlocks.end();
1224  // mit++
1225  // ) {
1226  // NonlinearElasticElement::CommonData &elastic_common_data = iso_elastic.commonData;
1227  // NonlinearElasticElement::BlockData &elastic_block_data = mit->second;
1228  // nitsche_element_iso.getOpPtrVector().push_back(
1229  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1230  // "DISPLACEMENT",elastic_common_data
1231  // )
1232  // );
1233  // nitsche_element_iso.getOpPtrVector().push_back(
1234  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1235  // "MESH_NODE_POSITIONS",elastic_common_data
1236  // )
1237  // );
1238  // nitsche_element_iso.getOpPtrVector().push_back(
1239  // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
1240  // "DISPLACEMENT",elastic_block_data,elastic_common_data,1,false,false,true
1241  // )
1242  // );
1243  // nitsche_element_iso.getOpPtrVector().push_back(
1244  // new PeriodicNitscheConstrains::OpCalculateHomogenisedStressSurfaceIntegral(
1245  // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
1246  // elastic_block_data,elastic_common_data,stress_homo
1247  // )
1248  // );
1249  // }
1250  // for(
1251  // map<int,NonlinearElasticElement::BlockData>::iterator mit = trans_elastic.setOfBlocks.begin();
1252  // mit!=trans_elastic.setOfBlocks.end();
1253  // mit++
1254  // ) {
1255  // NonlinearElasticElement::CommonData &elastic_common_data = trans_elastic.commonData;
1256  // NonlinearElasticElement::BlockData &elastic_block_data = mit->second;
1257  // nitsche_element_trans.getOpPtrVector().push_back(
1258  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1259  // "DISPLACEMENT",elastic_common_data
1260  // )
1261  // );
1262  // nitsche_element_trans.getOpPtrVector().push_back(
1263  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1264  // "MESH_NODE_POSITIONS",elastic_common_data
1265  // )
1266  // );
1267  // nitsche_element_trans.getOpPtrVector().push_back(
1268  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1269  // "POTENTIAL_FIELD",elastic_common_data
1270  // )
1271  // );
1272  // nitsche_element_trans.getOpPtrVector().push_back(
1273  // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
1274  // "DISPLACEMENT",elastic_block_data,elastic_common_data,2,false,false,true
1275  // )
1276  // );
1277  // nitsche_element_trans.getOpPtrVector().push_back(
1278  // new PeriodicNitscheConstrains::OpCalculateHomogenisedStressSurfaceIntegral(
1279  // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
1280  // elastic_block_data,elastic_common_data,stress_homo
1281  // )
1282  // );
1283  // }
1284  // } else {
1285  // for(
1286  // map<int,NonlinearElasticElement::BlockData>::iterator mit = iso_elastic.setOfBlocks.begin();
1287  // mit!=iso_elastic.setOfBlocks.end();
1288  // mit++
1289  // ) {
1290  // NonlinearElasticElement::CommonData &elastic_common_data = iso_elastic.commonData;
1291  // NonlinearElasticElement::BlockData &elastic_block_data = mit->second;
1292  // ave_stress_iso.getOpPtrVector().push_back(
1293  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1294  // "DISPLACEMENT",elastic_common_data
1295  // )
1296  // );
1297  // ave_stress_iso.getOpPtrVector().push_back(
1298  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1299  // "MESH_NODE_POSITIONS",elastic_common_data
1300  // )
1301  // );
1302  // ave_stress_iso.getOpPtrVector().push_back(
1303  // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
1304  // "DISPLACEMENT",elastic_block_data,elastic_common_data,1,false,false,true
1305  // )
1306  // );
1307  // ave_stress_iso.getOpPtrVector().push_back(
1308  // new PeriodicNitscheConstrains::OpCalculateHomogenisedStressVolumeIntegral(
1309  // "DISPLACEMENT",elastic_block_data,elastic_common_data,stress_homo
1310  // )
1311  // );
1312  // }
1313  // for(
1314  // map<int,NonlinearElasticElement::BlockData>::iterator mit = trans_elastic.setOfBlocks.begin();
1315  // mit!=trans_elastic.setOfBlocks.end();
1316  // mit++
1317  // ) {
1318  // NonlinearElasticElement::CommonData &elastic_common_data = trans_elastic.commonData;
1319  // NonlinearElasticElement::BlockData &elastic_block_data = mit->second;
1320  // ave_stress_trans.getOpPtrVector().push_back(
1321  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1322  // "DISPLACEMENT",elastic_common_data
1323  // )
1324  // );
1325  // ave_stress_trans.getOpPtrVector().push_back(
1326  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1327  // "MESH_NODE_POSITIONS",elastic_common_data
1328  // )
1329  // );
1330  // ave_stress_trans.getOpPtrVector().push_back(
1331  // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1332  // "POTENTIAL_FIELD",elastic_common_data
1333  // )
1334  // );
1335  // ave_stress_trans.getOpPtrVector().push_back(
1336  // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
1337  // "DISPLACEMENT",elastic_block_data,elastic_common_data,2,false,false,true
1338  // )
1339  // );
1340  // ave_stress_trans.getOpPtrVector().push_back(
1341  // new PeriodicNitscheConstrains::OpCalculateHomogenisedStressVolumeIntegral(
1342  // "DISPLACEMENT",elastic_block_data,elastic_common_data,stress_homo
1343  // )
1344  // );
1345  // }
1346  // }
1347  break;
1348  case NOHOMOBC:
1349  SETERRQ(PETSC_COMM_SELF,MOFEM_IMPOSIBLE_CASE,"*** Boundary conditions not set");
1350  }
1351 
1352  struct MyPostProc: public PostProcVolumeOnRefinedMesh {
1353 
1354  bool doPreProcess;
1355  bool doPostProcess;
1356 
1357  MyPostProc(MoFEM::Interface &m_field):
1358  PostProcVolumeOnRefinedMesh(m_field),
1359  doPreProcess(true),
1360  doPostProcess(true)
1361  {}
1362 
1363  void setDoPreProcess() { doPreProcess = true; }
1364  void unSetDoPreProcess() { doPreProcess = false; }
1365  void setDoPostProcess() { doPostProcess = true; }
1366  void unSetDoPostProcess() { doPostProcess = false; }
1367 
1368 
1369 
1370  PetscErrorCode preProcess() {
1371  PetscFunctionBegin;
1372  if(doPreProcess) {
1373  ierr = PostProcVolumeOnRefinedMesh::preProcess(); CHKERRQ(ierr);
1374  }
1375  PetscFunctionReturn(0);
1376  }
1377  PetscErrorCode postProcess() {
1378  PetscFunctionBegin;
1379  if(doPostProcess) {
1380  ierr = PostProcVolumeOnRefinedMesh::postProcess(); CHKERRQ(ierr);
1381  }
1382  PetscFunctionReturn(0);
1383  }
1384  };
1385 
1386  MyPostProc post_proc(m_field);
1387 
1388  ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
1389  ierr = post_proc.addFieldValuesPostProc("DISPLACEMENT"); CHKERRQ(ierr);
1390  ierr = post_proc.addFieldValuesGradientPostProc("DISPLACEMENT"); CHKERRQ(ierr);
1391  ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
1392  ierr = post_proc.addFieldValuesGradientPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
1393  if(trans_iso_blocks) {
1394  ierr = post_proc.addFieldValuesGradientPostProc("POTENTIAL_FIELD"); CHKERRQ(ierr);
1395  }
1396  for(
1397  map<int,NonlinearElasticElement::BlockData>::iterator sit = iso_elastic.setOfBlocks.begin();
1398  sit != iso_elastic.setOfBlocks.end(); sit++
1399  ) {
1400  post_proc.getOpPtrVector().push_back(
1401  new PostPorcStress(
1402  post_proc.postProcMesh,
1403  post_proc.mapGaussPts,
1404  "DISPLACEMENT",
1405  sit->second,
1406  post_proc.commonData,
1407  false
1408  )
1409  );
1410  }
1411  for(
1412  map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
1413  sit != trans_elastic.setOfBlocks.end(); sit++
1414  ) {
1415  post_proc.getOpPtrVector().push_back(
1416  new PostPorcStress(
1417  post_proc.postProcMesh,
1418  post_proc.mapGaussPts,
1419  "DISPLACEMENT",
1420  sit->second,
1421  post_proc.commonData
1422  )
1423  );
1424  }
1425 
1426  PetscScalar *avec;
1427  ierr = VecGetArray(stress_homo,&avec); CHKERRQ(ierr);
1428  for(int ii = 0;ii<6;ii++) {
1429  ierr = VecZeroEntries(D); CHKERRQ(ierr);
1430  ierr = KSPSolve(solver,F[ii],D); CHKERRQ(ierr);
1431  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1432  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1433  ierr = m_field.set_global_ghost_vector(
1434  "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
1435  ); CHKERRQ(ierr);
1436  post_proc.setDoPreProcess();
1437  post_proc.unSetDoPostProcess();
1438  ierr = m_field.loop_finite_elements(
1439  "ELASTIC_MECHANICS","ELASTIC",post_proc
1440  ); CHKERRQ(ierr);
1441  post_proc.unSetDoPreProcess();
1442  post_proc.setDoPostProcess();
1443 
1444  if(m_field.check_finite_element("TRAN_ISOTROPIC_ELASTIC"))
1445  ierr = m_field.loop_finite_elements(
1446  "ELASTIC_MECHANICS","TRAN_ISOTROPIC_ELASTIC",post_proc
1447  ); CHKERRQ(ierr);
1448  {
1449  ostringstream sss;
1450  sss << "mode_" << homo_bc_names[choise_value] << "_" << ii << ".h5m";
1451  rval = post_proc.postProcMesh.write_file(
1452  sss.str().c_str(),"MOAB","PARALLEL=WRITE_PART"
1453  ); CHKERRQ_MOAB(rval);
1454  }
1455  ierr = VecZeroEntries(stress_homo); CHKERRQ(ierr);
1456  if(choise_value == HOMOBCDISP) {
1457  ierr = m_field.loop_finite_elements(
1458  "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_disp.getLoopFeRVEBCStress()
1459  ); CHKERRQ(ierr);
1460  }
1461  if(choise_value == HOMOBCTRAC) {
1462  ierr = m_field.loop_finite_elements(
1463  "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_trac.getLoopFeRVEBCStress()
1464  ); CHKERRQ(ierr);
1465  }
1466  if(choise_value == HOMOBCPERIODIC) {
1467  ierr = m_field.loop_finite_elements(
1468  "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCStress()
1469  ); CHKERRQ(ierr);
1470  }
1471  // if(choise_value == NITSCHE) {
1472  // if(stress_by_boundary_integral) {
1473  // nitsche_element_iso.addToRule = 1;
1474  // nitsche_element_trans.addToRule = 1;
1475  // ierr = m_field.loop_finite_elements(
1476  // "ELASTIC_MECHANICS","ELASTIC",nitsche_element_iso
1477  // ); CHKERRQ(ierr);
1478  // if(m_field.check_finite_element("TRAN_ISOTROPIC_ELASTIC"))
1479  // ierr = m_field.loop_finite_elements(
1480  // "ELASTIC_MECHANICS","TRAN_ISOTROPIC_ELASTIC",nitsche_element_trans
1481  // ); CHKERRQ(ierr);
1482  // } else {
1483  // ave_stress_iso.addToRule = 1;
1484  // ave_stress_trans.addToRule = 1;
1485  // ierr = m_field.loop_finite_elements(
1486  // "ELASTIC_MECHANICS","ELASTIC",ave_stress_iso
1487  // ); CHKERRQ(ierr);
1488  // if(m_field.check_finite_element("TRAN_ISOTROPIC_ELASTIC"))
1489  // ierr = m_field.loop_finite_elements(
1490  // "ELASTIC_MECHANICS","TRAN_ISOTROPIC_ELASTIC",ave_stress_trans
1491  // ); CHKERRQ(ierr);
1492  // }
1493  // }
1494  ierr = PetscOptionsGetReal(
1495  PETSC_NULL,"-my_rve_volume",&rve_volume,PETSC_NULL
1496  ); CHKERRQ(ierr);
1497  ierr = VecAssemblyBegin(stress_homo); CHKERRQ(ierr);
1498  ierr = VecAssemblyEnd(stress_homo); CHKERRQ(ierr);
1499  ierr = VecScale(stress_homo,1.0/rve_volume); CHKERRQ(ierr);
1500  ierr = VecGhostUpdateBegin(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1501  ierr = VecGhostUpdateEnd(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1502  for(int jj=0; jj<6; jj++) {
1503  Dmat(jj,ii) = avec[jj];
1504  }
1505  }
1506  ierr = VecRestoreArray(stress_homo,&avec); CHKERRQ(ierr);
1507 
1508  PetscPrintf(
1509  PETSC_COMM_WORLD,"\nHomogenised Stiffens Matrix = \n\n"
1510  );
1511 
1512  for(int ii=0; ii<6; ii++) {
1513  PetscPrintf(
1514  PETSC_COMM_WORLD,
1515  "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",
1516  ii,Dmat(ii,0),Dmat(ii,1),Dmat(ii,2),Dmat(ii,3),Dmat(ii,4),Dmat(ii,5)
1517  );
1518  }
1519 
1520  //Saving Dmat as a bindary file to use it macro-structure
1521 
1522 
1523  char output_file_Dmat[255];
1524  ierr = PetscOptionsGetString(PETSC_NULL,"-my_output_file_Dmat",output_file_Dmat,255,&flg); CHKERRQ(ierr);
1525  if(flg) {
1526 
1527  //Reading and writing binary files
1528  if(pcomm->rank()==0){
1529  int fd;
1530  PetscViewer view_out;
1531  PetscViewerBinaryOpen(PETSC_COMM_SELF,output_file_Dmat,FILE_MODE_WRITE,&view_out);
1532  PetscViewerBinaryGetDescriptor(view_out,&fd);
1533  PetscBinaryWrite(fd,&Dmat(0,0),36,PETSC_DOUBLE,PETSC_FALSE);
1534  PetscViewerDestroy(&view_out);
1535  }
1536 
1537  // MatrixDouble Dmat1;
1538  // Dmat1.resize(6,6); Dmat1.clear();
1539  // if(pcomm->rank()==0){
1540  // int fd;
1541  // PetscViewer view_in;
1542  // PetscViewerBinaryOpen(PETSC_COMM_SELF,output_file_Dmat,FILE_MODE_READ,&view_in);
1543  // PetscViewerBinaryGetDescriptor(view_in,&fd);
1544  // PetscBinaryRead(fd,&Dmat1(0,0),36,PETSC_DOUBLE);
1545  // PetscViewerDestroy(&view_in);
1546  // cout<< "Dmat1 After Reading= "<<Dmat1<<endl;
1547  // }
1548  //
1549  }
1550 
1551  //detroy matrices
1552  for(int ii = 0;ii<6;ii++) {
1553  ierr = VecDestroy(&F[ii]); CHKERRQ(ierr);
1554  }
1555  ierr = VecDestroy(&D); CHKERRQ(ierr);
1556  ierr = MatDestroy(&Aij); CHKERRQ(ierr);
1557  ierr = KSPDestroy(&solver); CHKERRQ(ierr);
1558  ierr = VecDestroy(&stress_homo); CHKERRQ(ierr);
1559 
1560  ierr = PetscTime(&v2);CHKERRQ(ierr);
1561  ierr = PetscGetCPUTime(&t2);CHKERRQ(ierr);
1562 
1563  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Time = %f CPU Time = %f\n",pcomm->rank(),v2-v1,t2-t1);
1564 
1565  }
1566  CATCH_ERRORS;
1567 
1568 
1569  MoFEM::Core::Finalize();
1570 }
SIDESET
@ SIDESET
Definition: definitions.h:160
BCs_RVELagrange_Disp::getLoopFeRVEBCStress
MyTriFE & getLoopFeRVEBCStress()
Definition: BCs_RVELagrange_Disp.hpp:39
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
BCs_RVELagrange_Trac_Rigid_Rot
Definition: BCs_RVELagrange_Trac_Rigid_Rot.hpp:22
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.
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: rve_mechanical.cpp:134
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
BCs_RVELagrange_Periodic
Definition: BCs_RVELagrange_Periodic.hpp:19
MyPostProc
Definition: plot_base.cpp:34
BCs_RVELagrange_Disp.hpp
tetcircumcenter_tp
void tetcircumcenter_tp(double a[3], double b[3], double c[3], double d[3], double circumcenter[3], double *xi, double *eta, double *zeta)
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
NonlinearElasticElement::OpJacobianPiolaKirchhoffStress
Operator performs automatic differentiation.
Definition: NonLinearElasticElement.hpp:370
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:467
MoFEM::CoreInterface::modify_finite_element_add_field_row
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
NOBASE
@ NOBASE
Definition: definitions.h:59
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
MoFEM::Mat_Elastic
Elastic material data structure.
Definition: MaterialBlocks.hpp:139
main
int main(int argc, char *argv[])
Definition: rve_mechanical.cpp:136
NonlinearElasticElement::feRhs
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
Definition: NonLinearElasticElement.hpp:65
PostProcTemplateVolumeOnRefinedMesh< MoFEM::VolumeElementForcesAndSourcesCore >::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: PostProcOnRefMesh.hpp:797
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
BCs_RVELagrange_Disp::getLoopFeRVEBCRhs
MyTriFE & getLoopFeRVEBCRhs()
Definition: BCs_RVELagrange_Disp.hpp:36
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::get_field_meshset
virtual EntityHandle get_field_meshset(const std::string name) const =0
get field meshset
tricircumcenter3d_tp
void tricircumcenter3d_tp(double a[3], double b[3], double c[3], double circumcenter[3], double *xi, double *eta)
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:126
BCs_RVELagrange_Disp::setRVEBCsOperators
PetscErrorCode setRVEBCsOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &fvec)
Definition: BCs_RVELagrange_Disp.hpp:580
homo_bc_names
const char * homo_bc_names[]
Definition: rve_mechanical.cpp:87
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.
HOMOBCTRAC
@ HOMOBCTRAC
Definition: rve_mechanical.cpp:82
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
BCs_RVELagrange_Disp
Definition: BCs_RVELagrange_Disp.hpp:18
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
Composite_xyzcoord
Definition: homogenisation_periodic_atom_test.cpp:72
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
BCs_RVELagrange_Disp::setRVEBCsHomoStressOperators
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec Stress_Homo)
Definition: BCs_RVELagrange_Disp.hpp:671
NonlinearElasticElement::CommonData::spatialPositions
string spatialPositions
Definition: NonLinearElasticElement.hpp:109
xcoord_tag
Definition: homogenisation_periodic_atom_test.cpp:68
NonlinearElasticElement::feLhs
MyVolumeFE feLhs
Definition: NonLinearElasticElement.hpp:67
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
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
HOMOBCPERIODIC
@ HOMOBCPERIODIC
Definition: rve_mechanical.cpp:81
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::CoreInterface::modify_finite_element_add_field_col
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
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
a
constexpr double a
Definition: approx_sphere.cpp:30
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.
MoFEM::Mat_Elastic_TransIso::data
_data_ data
Definition: MaterialBlocks.hpp:382
Hooke.hpp
Implementation of linear elastic material.
MoFEM::Mat_Elastic_TransIso
Transverse Isotropic material data structure.
Definition: MaterialBlocks.hpp:371
BCs_RVELagrange_Trac_Rigid_Trans.hpp
MAT_ELASTICSET
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:172
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx
Definition: NonLinearElasticElement.hpp:556
COL
@ COL
Definition: definitions.h:136
HOMOBCDISP
@ HOMOBCDISP
Definition: rve_mechanical.cpp:80
eta
double eta
Definition: free_surface.cpp:170
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
HomoBCTypes
HomoBCTypes
Definition: rve_mechanical.cpp:79
NonlinearElasticElement::CommonData::meshPositions
string meshPositions
Definition: NonLinearElasticElement.hpp:110
NonLinearElasticElement.hpp
Operators and data structures for non-linear elastic analysis.
Projection10NodeCoordsOnField.hpp
BCs_RVELagrange_Trac
Definition: BCs_RVELagrange_Trac.hpp:22
MoFEM::Mat_Elastic::data
_data_ data
Definition: MaterialBlocks.hpp:155
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
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
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::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
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: rve_mechanical.cpp:102
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
BCs_RVELagrange_Trac::setRVEBCsOperators
PetscErrorCode setRVEBCsOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &f)
Definition: BCs_RVELagrange_Trac.hpp:383
NonlinearElasticElement::commonData
CommonData commonData
Definition: NonLinearElasticElement.hpp:121
NonlinearElasticElement::MyVolumeFE
definition of volume element
Definition: NonLinearElasticElement.hpp:31
Range
NOHOMOBC
@ NOHOMOBC
Definition: rve_mechanical.cpp:84
NITSCHE
@ NITSCHE
Definition: rve_mechanical.cpp:83
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
NonlinearElasticElement::OpGetCommonDataAtGaussPts
Definition: NonLinearElasticElement.hpp:362
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
NonlinearElasticElement::OpRhsPiolaKirchhoff
Definition: NonLinearElasticElement.hpp:520
MoFEM::DeprecatedCoreInterface::seed_ref_level
DEPRECATED MoFEMErrorCode seed_ref_level(const Range &ents, const BitRefLevel &bit, const bool only_tets=true, int verb=-1)
seed entities in the range and their adjacencies in a particular BitRefLevel
Definition: DeprecatedCoreInterface.cpp:24
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:955
PostPorcStress
DEPRECATED typedef PostProcStress PostPorcStress
Definition: PostProcStresses.hpp:193
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:71
BCs_RVELagrange_Trac_Rigid_Rot.hpp
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
BCs_RVELagrange_Periodic.hpp
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
SmallTransverselyIsotropic.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
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
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
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
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
BCs_RVELagrange_Trac_Rigid_Rot::setRVEBCsRigidBodyRotOperators
PetscErrorCode setRVEBCsRigidBodyRotOperators(string field_name, string lagrang_field_name, Mat aij, map< int, RVEBC_Data > setOfRVEBC)
Definition: BCs_RVELagrange_Trac_Rigid_Rot.hpp:105
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::check_finite_element
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
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
help
static char help[]
Definition: rve_mechanical.cpp:77
BCs_RVELagrange_Trac::setRVEBCsHomoStressOperators
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec stress_homo)
Definition: BCs_RVELagrange_Trac.hpp:448
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
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
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
zcoord_tag
Definition: homogenisation_periodic_atom_test.cpp:70
Generated by Doxygen 1.8.17 and hosted at University of Glasgow