v0.12.1
Classes | Typedefs | Enumerations | Functions | Variables
rve_mechanical.cpp File Reference

Calculates stiffness matrix for elastic RVE. More...

#include <MoFEM.hpp>
#include <Projection10NodeCoordsOnField.hpp>
#include <petsctime.h>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <MethodForForceScaling.hpp>
#include <TimeForceScale.hpp>
#include <BCs_RVELagrange_Disp.hpp>
#include <BCs_RVELagrange_Trac.hpp>
#include <BCs_RVELagrange_Trac_Rigid_Rot.hpp>
#include <BCs_RVELagrange_Trac_Rigid_Trans.hpp>
#include <BCs_RVELagrange_Periodic.hpp>
#include <adolc/adolc.h>
#include <NonLinearElasticElement.hpp>
#include <Hooke.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <SmallTransverselyIsotropic.hpp>
#include <VolumeCalculation.hpp>
#include <boost/ptr_container/ptr_map.hpp>
#include <PostProcOnRefMesh.hpp>
#include <PostProcStresses.hpp>
#include <moab/AdaptiveKDTree.hpp>
#include <algorithm>

Go to the source code of this file.

Classes

struct  Face_CenPos_Handle
 
struct  xcoord_tag
 
struct  ycoord_tag
 
struct  zcoord_tag
 
struct  Tri_Hand_tag
 
struct  Composite_xyzcoord
 

Typedefs

typedef 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
 

Enumerations

enum  HomoBCTypes {
  HOMOBCDISP , HOMOBCPERIODIC , HOMOBCTRAC , NITSCHE ,
  NOHOMOBC
}
 

Functions

void tetcircumcenter_tp (double a[3], double b[3], double c[3], double d[3], double circumcenter[3], double *xi, double *eta, double *zeta)
 
void tricircumcenter3d_tp (double a[3], double b[3], double c[3], double circumcenter[3], double *xi, double *eta)
 
int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 
const char * homo_bc_names []
 

Detailed Description

Calculates stiffness matrix for elastic RVE.

Three types of boundary conditions are implemented, i.e. HOMOBCDISP, HOMOBCPERIODIC, HOMOBCTRAC, NITSCHE.

NITSHCE method allow to apply periodic boundary conditions to arbitrary convex shape RVE.

Definition in file rve_mechanical.cpp.

Typedef Documentation

◆ Face_CenPos_Handle_multiIndex

typedef 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 at line 134 of file rve_mechanical.cpp.

Enumeration Type Documentation

◆ HomoBCTypes

enum HomoBCTypes
Enumerator
HOMOBCDISP 
HOMOBCPERIODIC 
HOMOBCTRAC 
NITSCHE 
NOHOMOBC 

Definition at line 79 of file rve_mechanical.cpp.

79  {
80  HOMOBCDISP,
81  HOMOBCPERIODIC,
82  HOMOBCTRAC,
83  NITSCHE,
84  NOHOMOBC
85 };
NOHOMOBC
@ NOHOMOBC
Definition: rve_mechanical.cpp:84
HOMOBCPERIODIC
@ HOMOBCPERIODIC
Definition: rve_mechanical.cpp:81
NITSCHE
@ NITSCHE
Definition: rve_mechanical.cpp:83
HOMOBCDISP
@ HOMOBCDISP
Definition: rve_mechanical.cpp:80
HOMOBCTRAC
@ HOMOBCTRAC
Definition: rve_mechanical.cpp:82

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 136 of file rve_mechanical.cpp.

136  {
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 }
PostPorcStress
DEPRECATED typedef PostProcStress PostPorcStress
Definition: PostProcStresses.hpp:206
COL
@ COL
Definition: definitions.h:136
ROW
@ ROW
Definition: definitions.h:136
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
NOBASE
@ NOBASE
Definition: definitions.h:72
MOAB_THROW
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:554
NOFIELD
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:97
H1
@ H1
continuous field
Definition: definitions.h:98
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MAT_ELASTICSET
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:172
SIDESET
@ SIDESET
Definition: definitions.h:160
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:467
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:46
MOFEM_IMPOSIBLE_CASE
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:48
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
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::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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
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_filed)=0
set finite element field data
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
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
MoFEM::DeprecatedCoreInterface::partition_ghost_dofs
DEPRECATED MoFEMErrorCode partition_ghost_dofs(const std::string &name, int verb=-1)
determine ghost nodes
Definition: DeprecatedCoreInterface.cpp:451
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.
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
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.
_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:61
_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:83
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::CoreInterface::add_problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
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
MoFEM::DeprecatedCoreInterface::get_entities_by_ref_level
DEPRECATED MoFEMErrorCode get_entities_by_ref_level(const BitRefLevel &bit, const BitRefLevel &mask, const EntityHandle meshset)
add all ents from ref level given by bit to meshset
Definition: DeprecatedCoreInterface.cpp:258
MoFEM::DeprecatedCoreInterface::get_entities_by_type_and_ref_level
DEPRECATED MoFEMErrorCode get_entities_by_type_and_ref_level(const BitRefLevel &bit, const BitRefLevel &mask, const EntityType type, const EntityHandle meshset, int verb=-1)
add all ents from ref level given by bit to meshset
Definition: DeprecatedCoreInterface.cpp:243
MoFEM::DeprecatedCoreInterface::set_global_ghost_vector
DEPRECATED MoFEMErrorCode set_global_ghost_vector(const Problem *problem_ptr, RowColData rc, Vec V, InsertMode mode, ScatterMode scatter_mode) const
set values of vector from/to mesh database (collective)
Definition: DeprecatedCoreInterface.cpp:196
MoFEM::DeprecatedCoreInterface::VecCreateGhost
DEPRECATED MoFEMErrorCode VecCreateGhost(const std::string &name, RowColData rc, Vec *V) const
create ghost vector for problem (collective)
Definition: DeprecatedCoreInterface.cpp:132
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
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:35
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:78
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:153
MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:163
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:193
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1095
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:214
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:183
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1949
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:32
help
static char help[]
Definition: rve_mechanical.cpp:77
homo_bc_names
const char * homo_bc_names[]
Definition: rve_mechanical.cpp:87
BCs_RVELagrange_Disp
Definition: BCs_RVELagrange_Disp.hpp:18
BCs_RVELagrange_Periodic
Definition: BCs_RVELagrange_Periodic.hpp:19
BCs_RVELagrange_Trac_Rigid_Rot
Definition: BCs_RVELagrange_Trac_Rigid_Rot.hpp:22
BCs_RVELagrange_Trac_Rigid_Trans
Definition: BCs_RVELagrange_Trac_Rigid_Trans.hpp:18
BCs_RVELagrange_Trac
Definition: BCs_RVELagrange_Trac.hpp:22
Composite_xyzcoord
Definition: homogenisation_periodic_atom_test.cpp:72
Face_CenPos_Handle
Definition: homogenisation_periodic_atom_test.cpp:61
Hooke
Hook equation.
Definition: Hooke.hpp:33
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
MoFEM::CoreInterface::get_field_meshset
virtual EntityHandle get_field_meshset(const std::string name) const =0
get field meshset
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::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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.
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:92
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:85
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:26
MoFEM::DeprecatedCoreInterface::MatCreateMPIAIJWithArrays
DEPRECATED MoFEMErrorCode MatCreateMPIAIJWithArrays(const std::string &name, Mat *Aij, int verb=DEFAULT_VERBOSITY)
Definition: DeprecatedCoreInterface.cpp:490
MoFEM::DeprecatedCoreInterface::get_cubit_msId_entities_by_dimension
DEPRECATED MoFEMErrorCode get_cubit_msId_entities_by_dimension(const int ms_id, const CubitBCType cubit_bc_type, const int dimension, Range &entities, const bool recursive=false)
Definition: DeprecatedCoreInterface.cpp:352
MoFEM::DeprecatedCoreInterface::seed_ref_level_MESHSET
DEPRECATED MoFEMErrorCode seed_ref_level_MESHSET(const EntityHandle meshset, const BitRefLevel &bit, int verb=-1)
Definition: DeprecatedCoreInterface.cpp:45
MoFEM::DeprecatedCoreInterface::seed_finite_elements
DEPRECATED MoFEMErrorCode seed_finite_elements(const EntityHandle meshset, int verb=-1)
Definition: DeprecatedCoreInterface.cpp:52
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:37
MoFEM::DeprecatedCoreInterface::partition_problem
DEPRECATED MoFEMErrorCode partition_problem(const std::string &name, int verb=-1)
partition problem dofs (collective)
Definition: DeprecatedCoreInterface.cpp:438
MoFEM::DeprecatedCoreInterface::partition_finite_elements
DEPRECATED MoFEMErrorCode partition_finite_elements(const std::string &name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=-1)
partition finite elements
Definition: DeprecatedCoreInterface.cpp:442
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:31
MoFEM::DeprecatedCoreInterface::loop_finite_elements
DEPRECATED MoFEMErrorCode loop_finite_elements(const Problem *problem_ptr, const std::string &fe_name, FEMethod &method, int lower_rank, int upper_rank, MoFEMTypes bh, int verb=DEFAULT_VERBOSITY)
Definition: DeprecatedCoreInterface.cpp:513
MoFEM::DeprecatedCoreInterface::build_problems
virtual DEPRECATED MoFEMErrorCode build_problems(int verb=-1)=0
build problem data structures
MoFEM::Mat_Elastic_TransIso
Transverse Isotropic material data structure.
Definition: MaterialBlocks.hpp:384
MoFEM::Mat_Elastic_TransIso::data
_data_ data
Definition: MaterialBlocks.hpp:395
MoFEM::Mat_Elastic
Elastic material data structure.
Definition: MaterialBlocks.hpp:152
MoFEM::Mat_Elastic::data
_data_ data
Definition: MaterialBlocks.hpp:168
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:36
MyPostProc
Definition: plot_base.cpp:52
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx
Definition: NonLinearElasticElement.hpp:570
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:40
PostProcTemplateVolumeOnRefinedMesh< MoFEM::VolumeElementForcesAndSourcesCore >::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: PostProcOnRefMesh.hpp:722
PostProcTemplateVolumeOnRefinedMesh< MoFEM::VolumeElementForcesAndSourcesCore >::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: PostProcOnRefMesh.hpp:816
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:977
Tri_Hand_tag
Definition: homogenisation_periodic_atom_test.cpp:71
VolumeCalculation
Calculate volume.
Definition: VolumeCalculation.hpp:31
xcoord_tag
Definition: homogenisation_periodic_atom_test.cpp:68
ycoord_tag
Definition: homogenisation_periodic_atom_test.cpp:69
zcoord_tag
Definition: homogenisation_periodic_atom_test.cpp:70

◆ 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 
)

◆ tricircumcenter3d_tp()

void tricircumcenter3d_tp ( double  a[3],
double  b[3],
double  c[3],
double  circumcenter[3],
double *  xi,
double *  eta 
)

Variable Documentation

◆ help

char help[] = "...\n\n"
static

Definition at line 77 of file rve_mechanical.cpp.

◆ homo_bc_names

const char* homo_bc_names[]
Initial value:
= {
"disp",
"periodic",
"trac",
"nitsche",
"nohomobc"
}

Definition at line 87 of file rve_mechanical.cpp.

Generated on Mon Jan 24 2022 21:34:48 for MoFEM by Doxygen 1.9.1 and hosted at University of Glasgow