v0.14.0
rve_mech_plas_periodic.cpp
Go to the documentation of this file.
1 /** \file small_strain_plasticity.cpp
2  * \ingroup small_strain_plasticity
3  * \brief Small strain plasticity example
4  *
5  */
6 
7 /* This file is part of MoFEM.
8  * MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20 
21 
22 #include <MoFEM.hpp>
23 using namespace MoFEM;
24 
25 #include <boost/program_options.hpp>
26 using namespace std;
27 namespace po = boost::program_options;
28 
29 #include <boost/numeric/ublas/vector_proxy.hpp>
30 #include <boost/numeric/ublas/matrix.hpp>
31 #include <boost/numeric/ublas/matrix_proxy.hpp>
32 #include <boost/numeric/ublas/vector.hpp>
33 #include <boost/numeric/ublas/symmetric.hpp>
34 
35 #include <MethodForForceScaling.hpp>
36 #include <TimeForceScale.hpp>
37 #include <VolumeCalculation.hpp>
38 #include <BCs_RVELagrange_Disp.hpp>
39 #include <BCs_RVELagrange_Trac.hpp>
40 #include "BCs_RVELagrange_Periodic.hpp"
41 #include <BCs_RVELagrange_Trac_Rigid_Trans.hpp>
42 
43 #ifndef WITH_MODULE_SMALL_STRAIN_PLASTICITY
44  #error "Install module small strain plasticity https://bitbucket.org/likask/mofem_um_small_strain_plasticity"
45 #endif
46 
47 #include <adolc/adolc.h>
48 #include <SmallStrainPlasticity.hpp>
49 #include <SmallStrainPlasticityMaterialModels.hpp>
50 
51 #include <Projection10NodeCoordsOnField.hpp>
52 #include <PostProcOnRefMesh.hpp>
53 
54 #include <string>
55 
56 using namespace boost::numeric;
57 
58 
59 
60 
61 static char help[] =
62 "...\n"
63 "\n";
64 
65 
66 
67 
68 //=================================================================================================================================
69 //Define class and multindex container to store data for traiangles on the boundary of the RVE (it cannot be defined within main)
70 //=================================================================================================================================
71 
72 struct Face_CenPos_Handle {
73  double xcoord, ycoord, zcoord;
74  const EntityHandle Tri_Hand;
75  Face_CenPos_Handle(double _xcoord, double _ycoord, double _zcoord, const EntityHandle _Tri_Hand):xcoord(_xcoord),
76  ycoord(_ycoord), zcoord(_zcoord), Tri_Hand(_Tri_Hand) {}
77 };
78 
79 struct xcoord_tag {}; //tags to used in the multindex container
80 struct ycoord_tag {};
81 struct zcoord_tag {};
82 struct Tri_Hand_tag {};
83 struct Composite_xyzcoord {};
84 
85 typedef multi_index_container<
86 Face_CenPos_Handle,
87 indexed_by<
88 ordered_non_unique<
89 tag<xcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord> >,
90 
91 ordered_non_unique<
92 tag<ycoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord> >,
93 
94 ordered_non_unique<
95 tag<zcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> >,
96 
97 ordered_unique<
98 tag<Tri_Hand_tag>, member<Face_CenPos_Handle,const EntityHandle,&Face_CenPos_Handle::Tri_Hand> >,
99 
100 ordered_unique<
101 tag<Composite_xyzcoord>,
102 composite_key<
103 Face_CenPos_Handle,
104 member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord>,
105 member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord>,
106 member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> > >
107 > > Face_CenPos_Handle_multiIndex;
108 
109 //============================================================================================
110 //============================================================================================
111 
112 
113 
114 
115 
116 int main(int argc, char *argv[]) {
117  MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
118 
119  try {
120 
121  moab::Core mb_instance;
122  moab::Interface& moab = mb_instance;
123  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
124  if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
125 
126  PetscBool flg = PETSC_TRUE;
127  char mesh_file_name[255];
128  ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
129  if(flg != PETSC_TRUE) {
130  SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
131  }
132 
133  PetscInt order;
134  ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
135  if(flg != PETSC_TRUE) {
136  order = 2;
137  }
138  PetscInt bubble_order;
139  ierr = PetscOptionsGetInt(PETSC_NULL,"-my_bubble_order",&bubble_order,&flg); CHKERRQ(ierr);
140  if(flg != PETSC_TRUE) {
141  bubble_order = order;
142  }
143 
144  // use this if your mesh is partitioned and you run code on parts,
145  // you can solve very big problems
146  PetscBool is_partitioned = PETSC_FALSE;
147  ierr = PetscOptionsGetBool(PETSC_NULL,"-my_is_partitioned",&is_partitioned,&flg); CHKERRQ(ierr);
148 
149  //Applied strain on the RVE (vector of length 6) strain=[xx, yy, zz, xy, xz, zy]^T
150  double mygiven_strain[6];
151  int nmax=6;
152  ierr = PetscOptionsGetRealArray(PETSC_NULL,"-my_given_strain",mygiven_strain,&nmax,&flg); CHKERRQ(ierr);
153  VectorDouble given_strain;
154  given_strain.resize(6);
155  cblas_dcopy(6, &mygiven_strain[0], 1, &given_strain(0), 1);
156  cout<<"given_strain ="<<given_strain<<endl;
157 
158  //Read mesh to MOAB
159  if(is_partitioned == PETSC_TRUE) {
160  const char *option;
161  option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION=PARALLEL_PARTITION;";
162  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
163  } else {
164  const char *option;
165  option = "";
166  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
167  }
168 
169 
170  //Create MoFEM (Joseph) database
171  MoFEM::Core core(moab);
172  MoFEM::Interface& m_field = core;
173 
174  //ref meshset ref level 0
175  ierr = m_field.seed_ref_level_3D(0,0); CHKERRQ(ierr);
176 
177  // stl::bitset see for more details
178  BitRefLevel bit_level0;
179  bit_level0.set(0);
180  EntityHandle meshset_level0;
181  rval = moab.create_meshset(MESHSET_SET,meshset_level0); CHKERRQ_MOAB(rval);
182  ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
183  ierr = m_field.get_entities_by_ref_level(bit_level0,BitRefLevel().set(),meshset_level0); CHKERRQ(ierr);
184 
185 
186 
187  //=======================================================================================================
188  //Add Periodic Prisims Between Triangles on -ve and +ve faces to implement periodic bounary conditions
189  //=======================================================================================================
190 
191  //Populating the Multi-index container with -ve triangles
192  Range SurTrisNeg;
193  ierr = m_field.get_cubit_msId_entities_by_dimension(101,SIDESET,2,SurTrisNeg,true); CHKERRQ(ierr);
194  ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
195  Face_CenPos_Handle_multiIndex Face_CenPos_Handle_varNeg, Face_CenPos_Handle_varPos;
196  double TriCen[3], coords_Tri[9];
197  for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
198  const EntityHandle* conn_face; int num_nodes_Tri;
199 
200  //get nodes attached to the face
201  rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
202  //get nodal coordinates
203  rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
204 
205  //Find triangle centriod
206  TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
207  TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
208  TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
209 
210  //round values to 3 disimal places
211  if(TriCen[0]>=0) TriCen[0]=double(int(TriCen[0]*1000.0+0.5))/1000.0; else TriCen[0]=double(int(TriCen[0]*1000.0-0.5))/1000.0; //-ve and +ve value
212  if(TriCen[1]>=0) TriCen[1]=double(int(TriCen[1]*1000.0+0.5))/1000.0; else TriCen[1]=double(int(TriCen[1]*1000.0-0.5))/1000.0;
213  if(TriCen[2]>=0) TriCen[2]=double(int(TriCen[2]*1000.0+0.5))/1000.0; else TriCen[2]=double(int(TriCen[2]*1000.0-0.5))/1000.0;
214 
215  // cout<<"\n\n\nTriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
216  //fill the multi-index container with centriod coordinates and triangle handles
217  Face_CenPos_Handle_varNeg.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
218  // for(int ii=0; ii<3; ii++) cout<<"TriCen "<<TriCen[ii]<<endl;
219  // cout<<endl<<endl;
220  }
221 
222  // double aaa;
223  // aaa=0.5011;
224  // cout<<"\n\n\n\n\nfloor(aaa+0.5) = "<<double(int(aaa*1000.0+0.5))/1000.0<<endl<<endl<<endl<<endl;
225  // aaa=-0.5011;
226  // cout<<"\n\n\n\n\nfloor(aaa+0.5) = "<<double(int(aaa*1000.0-0.5))/1000.0<<endl<<endl<<endl<<endl;
227 
228 
229  //Populating the Multi-index container with +ve triangles
230  Range SurTrisPos;
231  ierr = m_field.get_cubit_msId_entities_by_dimension(102,SIDESET,2,SurTrisPos,true); CHKERRQ(ierr);
232  ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 102 = %d\n",SurTrisPos.size()); CHKERRQ(ierr);
233  for(Range::iterator it = SurTrisPos.begin(); it!=SurTrisPos.end(); it++) {
234  const EntityHandle* conn_face; int num_nodes_Tri;
235 
236  //get nodes attached to the face
237  rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
238  //get nodal coordinates
239  rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
240 
241  //Find triangle centriod
242  TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
243  TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
244  TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
245 
246  //round values to 3 disimal places
247  if(TriCen[0]>=0) TriCen[0]=double(int(TriCen[0]*1000.0+0.5))/1000.0; else TriCen[0]=double(int(TriCen[0]*1000.0-0.5))/1000.0;
248  if(TriCen[1]>=0) TriCen[1]=double(int(TriCen[1]*1000.0+0.5))/1000.0; else TriCen[1]=double(int(TriCen[1]*1000.0-0.5))/1000.0;
249  if(TriCen[2]>=0) TriCen[2]=double(int(TriCen[2]*1000.0+0.5))/1000.0; else TriCen[2]=double(int(TriCen[2]*1000.0-0.5))/1000.0;
250  // cout<<"\n\n\nTriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
251 
252  //fill the multi-index container with centriod coordinates and triangle handles
253  Face_CenPos_Handle_varPos.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
254  }
255 
256  //Find minimum and maximum X, Y and Z coordinates of the RVE (using multi-index container)
257  double XcoordMin, YcoordMin, ZcoordMin, XcoordMax, YcoordMax, ZcoordMax;
258  typedef Face_CenPos_Handle_multiIndex::index<xcoord_tag>::type::iterator Tri_Xcoord_iterator;
259  typedef Face_CenPos_Handle_multiIndex::index<ycoord_tag>::type::iterator Tri_Ycoord_iterator;
260  typedef Face_CenPos_Handle_multiIndex::index<zcoord_tag>::type::iterator Tri_Zcoord_iterator;
261  Tri_Xcoord_iterator XcoordMin_it, XcoordMax_it;
262  Tri_Ycoord_iterator YcoordMin_it, YcoordMax_it;
263  Tri_Zcoord_iterator ZcoordMin_it, ZcoordMax_it;
264 
265  //XcoordMax_it-- because .end() will point iterator after the data range but .begin() will point the iteratore to the first value of range
266  XcoordMin_it=Face_CenPos_Handle_varNeg.get<xcoord_tag>().begin(); XcoordMin=XcoordMin_it->xcoord;
267  XcoordMax_it=Face_CenPos_Handle_varPos.get<xcoord_tag>().end(); XcoordMax_it--; XcoordMax=XcoordMax_it->xcoord;
268  YcoordMin_it=Face_CenPos_Handle_varNeg.get<ycoord_tag>().begin(); YcoordMin=YcoordMin_it->ycoord;
269  YcoordMax_it=Face_CenPos_Handle_varPos.get<ycoord_tag>().end(); YcoordMax_it--; YcoordMax=YcoordMax_it->ycoord;
270  ZcoordMin_it=Face_CenPos_Handle_varNeg.get<zcoord_tag>().begin(); ZcoordMin=ZcoordMin_it->zcoord;
271  ZcoordMax_it=Face_CenPos_Handle_varPos.get<zcoord_tag>().end(); ZcoordMax_it--; ZcoordMax=ZcoordMax_it->zcoord;
272 
273 // cout<<"XcoordMin "<<XcoordMin << " XcoordMax "<<XcoordMax <<endl;
274 // cout<<"YcoordMin "<<YcoordMin << " YcoordMax "<<YcoordMax <<endl;
275 // cout<<"ZcoordMin "<<ZcoordMin << " ZcoordMax "<<ZcoordMax <<endl;
276  double LxRVE, LyRVE, LzRVE;
277  LxRVE=XcoordMax-XcoordMin;
278  LyRVE=YcoordMax-YcoordMin;
279  LzRVE=ZcoordMax-ZcoordMin;
280 
281 
282  //Creating Prisims between triangles on -ve and +ve faces
283  typedef Face_CenPos_Handle_multiIndex::index<Tri_Hand_tag>::type::iterator Tri_Hand_iterator;
284  Tri_Hand_iterator Tri_Neg;
285  typedef Face_CenPos_Handle_multiIndex::index<Composite_xyzcoord>::type::iterator xyzcoord_iterator;
286  xyzcoord_iterator Tri_Pos;
287  Range PrismRange;
288  double XPos, YPos, ZPos;
289  //int count=0;
290 
291  //loop over -ve triangles to create prisims elemenet between +ve and -ve triangles
292  for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
293 
294  Tri_Neg=Face_CenPos_Handle_varNeg.get<Tri_Hand_tag>().find(*it);
295  //cout<<"xcoord= "<<Tri_iit->xcoord << " ycoord= "<< Tri_iit->ycoord << " ycoord= "<< Tri_iit->zcoord <<endl;
296 
297  //corresponding +ve triangle
298  if(Tri_Neg->xcoord==XcoordMin){XPos=XcoordMax; YPos=Tri_Neg->ycoord; ZPos=Tri_Neg->zcoord;};
299  if(Tri_Neg->ycoord==YcoordMin){XPos=YPos=Tri_Neg->xcoord; YPos=YcoordMax; ZPos=Tri_Neg->zcoord;};
300  if(Tri_Neg->zcoord==ZcoordMin){XPos=YPos=Tri_Neg->xcoord; YPos=Tri_Neg->ycoord; ZPos=ZcoordMax; };
301  Tri_Pos=Face_CenPos_Handle_varPos.get<Composite_xyzcoord>().find(boost::make_tuple(XPos, YPos, ZPos));
302  // cout<<"Tri_Neg->xcoord= "<<Tri_Neg->xcoord << " Tri_Neg->ycoord "<< Tri_Neg->ycoord << " Tri_Neg->zcoord= "<< Tri_Neg->zcoord <<endl;
303  // cout<<"Tri_Pos->xcoord= "<<Tri_Pos->xcoord << " Tri_Pos->ycoord "<< Tri_Pos->ycoord << " Tri_Pos->zcoord= "<< Tri_Pos->zcoord <<endl;
304 
305 
306  //+ve and -ve nodes and their coords (+ve and -ve tiangles nodes can have matching problems, which can produce twisted prism)
307  EntityHandle PrismNodes[6];
308  vector<EntityHandle> TriNodesNeg, TriNodesPos;
309  double CoordNodeNeg[9], CoordNodePos[9];
310  rval = moab.get_connectivity(&(Tri_Neg->Tri_Hand),1,TriNodesNeg,true); CHKERRQ_MOAB(rval);
311  rval = moab.get_connectivity(&(Tri_Pos->Tri_Hand),1,TriNodesPos,true); CHKERRQ_MOAB(rval);
312  rval = moab.get_coords(&TriNodesNeg[0],3,CoordNodeNeg); MOAB_THROW(rval);
313  rval = moab.get_coords(&TriNodesPos[0],3,CoordNodePos); MOAB_THROW(rval);
314  for(int ii=0; ii<3; ii++){
315  PrismNodes[ii]=TriNodesNeg[ii];
316  }
317  // for(int ii=0; ii<3; ii++){
318  // cout<<"xcoord= "<<CoordNodeNeg[3*ii] << " ycoord= "<< CoordNodeNeg[3*ii+1] << " zcoord= "<< CoordNodeNeg[3*ii+2] <<endl;
319  // }
320  // for(int ii=0; ii<3; ii++){
321  // cout<<"xcoord= "<<CoordNodePos[3*ii] << " ycoord= "<< CoordNodePos[3*ii+1] << " zcoord= "<< CoordNodePos[3*ii+2] <<endl;
322  // }
323 
324  //Match exact nodes to each other to avoide the problem of twisted prisms
325  double XNodeNeg, YNodeNeg, ZNodeNeg, XNodePos, YNodePos, ZNodePos;
326  for(int ii=0; ii<3; ii++){
327  if(Tri_Neg->xcoord==XcoordMin){XNodeNeg=XcoordMax; YNodeNeg=CoordNodeNeg[3*ii+1]; ZNodeNeg=CoordNodeNeg[3*ii+2];};
328  if(Tri_Neg->ycoord==YcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=YcoordMax; ZNodeNeg=CoordNodeNeg[3*ii+2];};
329  if(Tri_Neg->zcoord==ZcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=CoordNodeNeg[3*ii+1]; ZNodeNeg=ZcoordMax;};
330  for(int jj=0; jj<3; jj++){
331  //Round nodal coordinates to 3 dicimal places only for comparison
332  //round values to 3 disimal places
333  if(XNodeNeg>=0) XNodeNeg=double(int(XNodeNeg*1000.0+0.5))/1000.0; else XNodeNeg=double(int(XNodeNeg*1000.0-0.5))/1000.0;
334  if(YNodeNeg>=0) YNodeNeg=double(int(YNodeNeg*1000.0+0.5))/1000.0; else YNodeNeg=double(int(YNodeNeg*1000.0-0.5))/1000.0;
335  if(ZNodeNeg>=0) ZNodeNeg=double(int(ZNodeNeg*1000.0+0.5))/1000.0; else ZNodeNeg=double(int(ZNodeNeg*1000.0-0.5))/1000.0;
336 
337  XNodePos=CoordNodePos[3*jj]; YNodePos=CoordNodePos[3*jj+1]; ZNodePos=CoordNodePos[3*jj+2];
338  if(XNodePos>=0) XNodePos=double(int(XNodePos*1000.0+0.5))/1000.0; else XNodePos=double(int(XNodePos*1000.0-0.5))/1000.0;
339  if(YNodePos>=0) YNodePos=double(int(YNodePos*1000.0+0.5))/1000.0; else YNodePos=double(int(YNodePos*1000.0-0.5))/1000.0;
340  if(ZNodePos>=0) ZNodePos=double(int(ZNodePos*1000.0+0.5))/1000.0; else ZNodePos=double(int(ZNodePos*1000.0-0.5))/1000.0;
341 
342  if(XNodeNeg==XNodePos && YNodeNeg==YNodePos && ZNodeNeg==ZNodePos){
343  PrismNodes[3+ii]=TriNodesPos[jj];
344  break;
345  }
346  }
347  }
348  //prism nodes and their coordinates
349  double CoordNodesPrisms[18];
350  rval = moab.get_coords(PrismNodes,6,CoordNodesPrisms); MOAB_THROW(rval);
351  // for(int ii=0; ii<6; ii++){
352  // cout<<"xcoord= "<<CoordNodesPrisms[3*ii] << " ycoord= "<< CoordNodesPrisms[3*ii+1] << " zcoord= "<< CoordNodesPrisms[3*ii+2] <<endl;
353  // }
354  // cout<<endl<<endl;
355  //insertion of individula prism element and its addition to range PrismRange
356  EntityHandle PeriodicPrism;
357  rval = moab.create_element(MBPRISM,PrismNodes,6,PeriodicPrism); CHKERRQ_MOAB(rval);
358  PrismRange.insert(PeriodicPrism);
359 
360  //Adding Prisims to Element LAGRANGE_ELEM (to loop over these prisims)
361  EntityHandle PrismRangeMeshset;
362  rval = moab.create_meshset(MESHSET_SET,PrismRangeMeshset); CHKERRQ_MOAB(rval);
363  rval = moab.add_entities(PrismRangeMeshset,PrismRange); CHKERRQ_MOAB(rval);
364  ierr = m_field.seed_ref_level_3D(PrismRangeMeshset,bit_level0); CHKERRQ(ierr);
365  ierr = m_field.get_entities_by_ref_level(bit_level0,BitRefLevel().set(),meshset_level0); CHKERRQ(ierr);
366 
367  // //to see individual prisms
368  // Range Prism1;
369  // Prism1.insert(PeriodicPrism);
370  // EntityHandle out_meshset1;
371  // rval = moab.create_meshset(MESHSET_SET,out_meshset1); CHKERRQ_MOAB(rval);
372  // rval = moab.add_entities(out_meshset1,Prism1); CHKERRQ_MOAB(rval);
373  // ostringstream sss;
374  // sss << "Prism" << count << ".vtk"; count++;
375  // rval = moab.write_file(sss.str().c_str(),"VTK","",&out_meshset1,1); CHKERRQ_MOAB(rval);
376 
377  }
378 
379 
380  //cout<<"PrismRange "<<PrismRange<<endl;
381  // //Saving prisms in interface.vtk
382  // EntityHandle out_meshset1;
383  // rval = moab.create_meshset(MESHSET_SET,out_meshset1); CHKERRQ_MOAB(rval);
384  // rval = moab.add_entities(out_meshset1,PrismRange); CHKERRQ_MOAB(rval);
385  // rval = moab.write_file("Prisms.vtk","VTK","",&out_meshset1,1); CHKERRQ_MOAB(rval);
386 
387 
388 // //=======================================================================================================
389 // //=======================================================================================================
390 
391 
392 
393  SmallStrainParaboloidalPlasticity cp;
394  {
395  cp.tAgs.resize(3);
396  cp.tAgs[0] = 0;
397  cp.tAgs[1] = 1;
398  cp.tAgs[2] = 2;
399  cp.tOl = 1e-12;
400 
401  double young_modulus = 1;
402  double poisson_ratio = 0.25;
403  cp.sIgma_yt = 1;
404  cp.sIgma_yc = 1;
405 
406  cp.Ht = 0.1;
407  cp.Hc = 0.1;
408 
409  cp.nup = 0.3;
410 
411 
412  {
413  ierr = PetscOptionsGetReal(0,"-my_young_modulus",&young_modulus,0); CHKERRQ(ierr);
414  ierr = PetscOptionsGetReal(0,"-my_poisson_ratio",&poisson_ratio,0); CHKERRQ(ierr);
415  cp.mu = MU(young_modulus,poisson_ratio);
416  cp.lambda = LAMBDA(young_modulus,poisson_ratio);
417  ierr = PetscOptionsGetReal(0,"-my_sigma_yt",&cp.sIgma_yt,0); CHKERRQ(ierr);
418  ierr = PetscOptionsGetReal(0,"-my_sigma_yc",&cp.sIgma_yc,0); CHKERRQ(ierr);
419 
420  ierr = PetscOptionsGetReal(0,"-my_Ht",&cp.Ht,0); CHKERRQ(ierr);
421  ierr = PetscOptionsGetReal(0,"-my_Hc",&cp.Hc,0); CHKERRQ(ierr);
422 
423  ierr = PetscOptionsGetReal(0,"-my_nt",&cp.nt,0); CHKERRQ(ierr);
424  ierr = PetscOptionsGetReal(0,"-my_nc",&cp.nc,0); CHKERRQ(ierr);
425 
426 
427  ierr = PetscOptionsGetReal(0,"-my_nup",&cp.nup,0); CHKERRQ(ierr);
428  }
429 
430  PetscPrintf(PETSC_COMM_WORLD,"young_modulus = %4.2e \n",young_modulus);
431  PetscPrintf(PETSC_COMM_WORLD,"poisson_ratio = %4.2e \n",poisson_ratio);
432 
433  PetscPrintf(PETSC_COMM_WORLD,"sIgma_yt = %4.2e \n",cp.sIgma_yt);
434  PetscPrintf(PETSC_COMM_WORLD,"sIgma_yc = %4.2e \n",cp.sIgma_yt);
435 
436  PetscPrintf(PETSC_COMM_WORLD,"nup = %4.2e \n",cp.nup);
437 
438  cp.sTrain.resize(6,false);
439  cp.sTrain.clear();
440  cp.plasticStrain.resize(6,false);
441  cp.plasticStrain.clear();
442  cp.internalVariables.resize(2,false);
443  cp.internalVariables.clear();
444  cp.createMatAVecR();
445  cp.snesCreate();
446  // ierr = SNESSetFromOptions(cp.sNes); CHKERRQ(ierr);
447  }
448 
449 
450 
451 // SmallStrainJ2Plasticity cp;
452 // {
453 // cp.tAgs.resize(3);
454 // cp.tAgs[0] = 0;
455 // cp.tAgs[1] = 1;
456 // cp.tAgs[2] = 2;
457 // cp.tOl = 1e-12;
458 //
459 // double young_modulus = 1;
460 // double poisson_ratio = 0.25;
461 // cp.sIgma_y = 1;
462 // cp.H = 0.1;
463 // cp.K = 0;
464 // cp.phi = 1;
465 // {
466 // ierr = PetscOptionsGetReal(0,"-my_young_modulus",&young_modulus,0); CHKERRQ(ierr);
467 // ierr = PetscOptionsGetReal(0,"-my_poisson_ratio",&poisson_ratio,0); CHKERRQ(ierr);
468 // cp.mu = MU(young_modulus,poisson_ratio);
469 // cp.lambda = LAMBDA(young_modulus,poisson_ratio);
470 // ierr = PetscOptionsGetReal(0,"-my_sigma_y",&cp.sIgma_y,0); CHKERRQ(ierr);
471 // ierr = PetscOptionsGetReal(0,"-my_H",&cp.H,0); CHKERRQ(ierr);
472 // ierr = PetscOptionsGetReal(0,"-my_K",&cp.K,0); CHKERRQ(ierr);
473 // ierr = PetscOptionsGetReal(0,"-my_phi",&cp.phi,0); CHKERRQ(ierr);
474 // }
475 //
476 // PetscPrintf(PETSC_COMM_WORLD,"young_modulus = %4.2e \n",young_modulus);
477 // PetscPrintf(PETSC_COMM_WORLD,"poisson_ratio = %4.2e \n",poisson_ratio);
478 // PetscPrintf(PETSC_COMM_WORLD,"sIgma_y = %4.2e \n",cp.sIgma_y);
479 // PetscPrintf(PETSC_COMM_WORLD,"H = %4.2e \n",cp.H);
480 // PetscPrintf(PETSC_COMM_WORLD,"K = %4.2e \n",cp.K);
481 // PetscPrintf(PETSC_COMM_WORLD,"phi = %4.2e \n",cp.phi);
482 //
483 //
484 // cp.sTrain.resize(6,false);
485 // cp.sTrain.clear();
486 // cp.plasticStrain.resize(6,false);
487 // cp.plasticStrain.clear();
488 // cp.internalVariables.resize(7,false);
489 // cp.internalVariables.clear();
490 // cp.createMatAVecR();
491 // cp.snesCreate();
492 // // ierr = SNESSetFromOptions(cp.sNes); CHKERRQ(ierr);
493 //
494 // }
495 
496 
497 // BitRefLevel problem_bit_level = bit_levels.back();
498 // Range CubitSideSets_meshsets;
499 // ierr = m_field.get_cubit_meshsets(SIDESET,CubitSideSets_meshsets); CHKERRQ(ierr);
500 
501  //Fields
502  ierr = m_field.add_field("DISPLACEMENT",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
503  ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
504  ierr = m_field.add_field("LAGRANGE_MUL_PERIODIC",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr); //For lagrange multipliers to control the periodic motion
505  ierr = m_field.add_field("LAGRANGE_MUL_RIGID_TRANS",NOFIELD,NOBASE,3); CHKERRQ(ierr); //To control the rigid body motion (3 Traslations and 3
506 
507 
508  //add entitities (by tets) to the field
509  ierr = m_field.add_ents_to_field_by_type(0,MBTET,"DISPLACEMENT"); CHKERRQ(ierr);
510  ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
511 
512 
513  Range surface_negative;
514  ierr = m_field.get_cubit_msId_entities_by_dimension(101,SIDESET,2,surface_negative,true); CHKERRQ(ierr);
515  ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 101 = %d\n",surface_negative.size()); CHKERRQ(ierr);
516  ierr = m_field.add_ents_to_field_by_type(surface_negative,MBTRI,"LAGRANGE_MUL_PERIODIC"); CHKERRQ(ierr);
517 
518 
519  //set app. order
520  ierr = m_field.set_field_order(0,MBTET,"DISPLACEMENT",bubble_order); CHKERRQ(ierr);
521  ierr = m_field.set_field_order(0,MBTRI,"DISPLACEMENT",order); CHKERRQ(ierr);
522  ierr = m_field.set_field_order(0,MBEDGE,"DISPLACEMENT",order); CHKERRQ(ierr);
523  ierr = m_field.set_field_order(0,MBVERTEX,"DISPLACEMENT",1); CHKERRQ(ierr);
524 
525  ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
526  ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
527  ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
528  ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
529 
530  ierr = m_field.set_field_order(0,MBTRI,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
531  ierr = m_field.set_field_order(0,MBEDGE,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
532  ierr = m_field.set_field_order(0,MBVERTEX,"LAGRANGE_MUL_PERIODIC",1); CHKERRQ(ierr);
533 
534 
535  //build field
536  ierr = m_field.build_fields(); CHKERRQ(ierr);
537  //10 node tets
538  Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
539  ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material,0); CHKERRQ(ierr);
540 
541  //FE
542  ierr = m_field.add_finite_element("PLASTIC"); CHKERRQ(ierr);
543  //Define rows/cols and element data
544  ierr = m_field.modify_finite_element_add_field_row("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
545  ierr = m_field.modify_finite_element_add_field_col("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
546  ierr = m_field.modify_finite_element_add_field_data("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
547  ierr = m_field.modify_finite_element_add_field_data("PLASTIC","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
548  ierr = m_field.add_ents_to_finite_element_by_type(0,MBTET,"PLASTIC"); CHKERRQ(ierr);
549 
550  BCs_RVELagrange_Periodic lagrangian_element_periodic(m_field);
551  BCs_RVELagrange_Trac lagrangian_element_trac(m_field);
552 
553  lagrangian_element_periodic.addLagrangiangElement("LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",PrismRange);
554  lagrangian_element_trac.addLagrangiangElement("LAGRANGE_ELEM_TRANS","DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS","MESH_NODE_POSITIONS");
555 
556  //build finite elements
557  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
558  //build adjacencies
559  ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
560 
561 
562  DMType dm_name = "PLASTIC_PROB";
563  ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
564 
565  DM dm;
566  ierr = DMCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(ierr);
567  ierr = DMSetType(dm,dm_name);CHKERRQ(ierr);
568 
569  //set dm datastruture which created mofem datastructures
570  ierr = DMMoFEMCreateMoFEM(dm,&m_field,dm_name,bit_level0); CHKERRQ(ierr);
571  ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
572  ierr = DMMoFEMSetIsPartitioned(dm,is_partitioned); CHKERRQ(ierr);
573  //add elements to dm
574  ierr = DMMoFEMAddElement(dm,"PLASTIC"); CHKERRQ(ierr);
575  ierr = DMMoFEMAddElement(dm,"LAGRANGE_ELEM"); CHKERRQ(ierr);
576  ierr = DMMoFEMAddElement(dm,"LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
577  ierr = DMSetUp(dm); CHKERRQ(ierr);
578 
579  //create matrices
580  Vec F,D;
581  Mat Aij;
582  ierr = DMCreateGlobalVector_MoFEM(dm,&D); CHKERRQ(ierr);
583  ierr = VecDuplicate(D,&F); CHKERRQ(ierr);
584  ierr = DMCreateMatrix_MoFEM(dm,&Aij); CHKERRQ(ierr);
585 
586  ierr = VecZeroEntries(D); CHKERRQ(ierr);
587  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
588  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
589  ierr = DMoFEMMeshToLocalVector(dm,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
590 
591  ierr = VecZeroEntries(F); CHKERRQ(ierr);
592  ierr = VecGhostUpdateBegin(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
593  ierr = VecGhostUpdateEnd(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
594  ierr = DMoFEMMeshToLocalVector(dm,F,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
595 
596  ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
597 
598  vector<Vec> Fvec(6); //jthis will be used to caluclate the homogenised stiffness matix
599  for(int ii = 0;ii<6;ii++) {
600  ierr = VecDuplicate(D,&Fvec[ii]); CHKERRQ(ierr);
601  ierr = VecZeroEntries(Fvec[ii]); CHKERRQ(ierr);
602  ierr = VecGhostUpdateBegin(Fvec[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
603  ierr = VecGhostUpdateEnd(Fvec[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
604  }
605 
606  //assemble Aij and F
607  SmallStrainPlasticity small_strain_plasticity(m_field);
608  {
609  PetscBool bbar = PETSC_TRUE;
610  ierr = PetscOptionsGetBool(0,"-my_bbar",&bbar,0); CHKERRQ(ierr);
611  small_strain_plasticity.commonData.bBar = bbar;
612  }
613  {
614  small_strain_plasticity.feRhs.getOpPtrVector().push_back(
615  new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
616  "DISPLACEMENT",small_strain_plasticity.commonData
617  )
618  );
619  small_strain_plasticity.feRhs.getOpPtrVector().push_back(
620  new SmallStrainPlasticity::OpCalculateStress(
621  m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp,false
622  )
623  );
624  small_strain_plasticity.feRhs.getOpPtrVector().push_back(
625  new SmallStrainPlasticity::OpAssembleRhs(
626  "DISPLACEMENT",small_strain_plasticity.commonData
627  )
628  );
629  small_strain_plasticity.feLhs.getOpPtrVector().push_back(
630  new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
631  "DISPLACEMENT",small_strain_plasticity.commonData
632  )
633  );
634  small_strain_plasticity.feLhs.getOpPtrVector().push_back(
635  new SmallStrainPlasticity::OpCalculateStress(
636  m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp,false
637  )
638  );
639  small_strain_plasticity.feLhs.getOpPtrVector().push_back(
640  new SmallStrainPlasticity::OpAssembleLhs(
641  "DISPLACEMENT",small_strain_plasticity.commonData
642  )
643  );
644  small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
645  new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
646  "DISPLACEMENT",small_strain_plasticity.commonData
647  )
648  );
649  small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
650  new SmallStrainPlasticity::OpCalculateStress(
651  m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp,false
652  )
653  );
654  small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
655  new SmallStrainPlasticity::OpUpdate(
656  "DISPLACEMENT",small_strain_plasticity.commonData
657  )
658  );
659  ierr = small_strain_plasticity.createInternalVariablesTag(); CHKERRQ(ierr);
660  }
661 
662  lagrangian_element_periodic.setRVEBCsOperatorsNonlinear("DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",Aij,Fvec,F,given_strain);
663 
664  BCs_RVELagrange_Trac_Rigid_Trans lagrangian_element_rigid_body_trans(m_field);
665  lagrangian_element_rigid_body_trans.setRVEBCsRigidBodyTranOperators("DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS",Aij, lagrangian_element_periodic.setOfRVEBC);
666 
667 
668  TimeForceScale time_force_scale("-my_macro_strian_history",false);
669  lagrangian_element_periodic.methodsOp.push_back(new TimeForceScale("-my_macro_strian_history",false));
670 
671 
672 
673 
674  //Adding elements to DMSnes
675  //Rhs
676  ierr = DMMoFEMSNESSetFunction(dm,"PLASTIC",&small_strain_plasticity.feRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
677  ierr = DMMoFEMSNESSetFunction(dm,"LAGRANGE_ELEM",&lagrangian_element_periodic.feRVEBCRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
678  ierr = DMMoFEMSNESSetFunction(dm,"LAGRANGE_ELEM",&lagrangian_element_periodic.feRVEBCRhsResidual,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
679 
680  //Lhs
681  ierr = DMMoFEMSNESSetJacobian(dm,"PLASTIC",&small_strain_plasticity.feLhs,NULL,NULL); CHKERRQ(ierr);
682  ierr = DMMoFEMSNESSetJacobian(dm,"LAGRANGE_ELEM",&lagrangian_element_periodic.feRVEBCLhs,NULL,NULL); CHKERRQ(ierr);
683  ierr = DMMoFEMSNESSetJacobian(dm,"LAGRANGE_ELEM_TRANS",&lagrangian_element_rigid_body_trans.feRVEBCLhs,NULL,NULL); CHKERRQ(ierr);
684 
685 
686  // Create Time Stepping solver
687  SNES snes;
688  SnesCtx *snes_ctx;
689  ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
690  //ierr = SNESSetDM(snes,dm); CHKERRQ(ierr);
691  ierr = DMMoFEMGetSnesCtx(dm,&snes_ctx); CHKERRQ(ierr);
692  ierr = SNESSetFunction(snes,F,SnesRhs,snes_ctx); CHKERRQ(ierr);
693  ierr = SNESSetJacobian(snes,Aij,Aij,SnesMat,snes_ctx); CHKERRQ(ierr);
694  ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
695 
696  PostProcVolumeOnRefinedMesh post_proc(m_field);
697  ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
698  ierr = post_proc.addFieldValuesPostProc("DISPLACEMENT"); CHKERRQ(ierr);
699  ierr = post_proc.addFieldValuesGradientPostProc("DISPLACEMENT"); CHKERRQ(ierr);
700  ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
701 
702  // Volume calculation
703  //==============================================================================================================================
704  VolumeElementForcesAndSourcesCore vol_elem(m_field);
705  Vec volume_vec;
706  int volume_vec_ghost[] = { 0 };
707  ierr = VecCreateGhost(
708  PETSC_COMM_WORLD,(!m_field.get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
709  ); CHKERRQ(ierr);
710  ierr = VecZeroEntries(volume_vec); CHKERRQ(ierr);
711  vol_elem.getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
712 
713  ierr = m_field.loop_finite_elements("PLASTIC_PROB","PLASTIC", vol_elem); CHKERRQ(ierr);
714 
715  ierr = VecAssemblyBegin(volume_vec); CHKERRQ(ierr);
716  ierr = VecAssemblyEnd(volume_vec); CHKERRQ(ierr);
717  double rve_volume;
718  ierr = VecSum(volume_vec,&rve_volume); CHKERRQ(ierr);
719 
720  // ierr = VecView(volume_vec,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
721  ierr = PetscPrintf(PETSC_COMM_WORLD,"RVE Volume %3.2g\n",rve_volume); CHKERRQ(ierr);
722  ierr = VecDestroy(&volume_vec);
723  //=============================================================================================================================
724 
725  double final_time = 1,delta_time = 0.1;
726  ierr = PetscOptionsGetReal(0,"-my_final_time",&final_time,0); CHKERRQ(ierr);
727  ierr = PetscOptionsGetReal(0,"-my_delta_time",&delta_time,0); CHKERRQ(ierr);
728  double delta_time0 = delta_time;
729 
730  // ierr = MatView(Aij,PETSC_VIEWER_DRAW_SELF); CHKERRQ(ierr);
731  // std::string wait;
732  // std::cin >> wait;
733 
734  Vec D0;
735  ierr = VecDuplicate(D,&D0); CHKERRQ(ierr);
736 
737  int step = 0;
738  double t = 0;
739  SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
740  for(;t<final_time;step++) {
741  t += delta_time;
742 
743  if(t>final_time){
744  break;
745  }
746  PetscPrintf(
747  PETSC_COMM_WORLD,"Step %d Time %6.4g final time %3.2g\n",step,t,final_time
748  );
749  //set time
750  lagrangian_element_periodic.getLoopFeRVEBCRhs().ts_t = t;
751  //solve problem at step
752  ierr = VecAssemblyBegin(D);
753  ierr = VecAssemblyEnd(D);
754  ierr = VecCopy(D,D0); CHKERRQ(ierr);
755 
756  if(step == 0 || reason < 0) {
757  ierr = SNESSetLagJacobian(snes,-2); CHKERRQ(ierr);
758  } else {
759  ierr = SNESSetLagJacobian(snes,-1); CHKERRQ(ierr);
760  }
761 
762 
763  ierr = SNESSolve(snes,PETSC_NULL,D); CHKERRQ(ierr);
764  int its;
765  ierr = SNESGetIterationNumber(snes,&its); CHKERRQ(ierr);
766 
767  ierr = PetscPrintf(
768  PETSC_COMM_WORLD,"number of Newton iterations = %D\n",its
769  ); CHKERRQ(ierr);
770 
771  ierr = SNESGetConvergedReason(snes,&reason); CHKERRQ(ierr);
772 // cout<<"reason ====== "<<reason<<endl;
773 
774  if(reason<0) {
775  t -= delta_time;
776  delta_time *= 0.1;
777  ierr = VecCopy(D0,D); CHKERRQ(ierr);
778  } else {const int its_d = 20;
779  const double gamma = 0.5;
780  const double max_reudction = 1;
781  const double min_reduction = 1e-1;
782  double reduction;
783  reduction = pow((double)its_d/(double)(its+1),gamma);
784  if(delta_time >= max_reudction*delta_time0 && reduction > 1) {
785  reduction = 1;
786  } else if(delta_time <= min_reduction*delta_time0 && reduction < 1) {
787  reduction = 1;
788  }
789 
790  ierr = PetscPrintf(
791  PETSC_COMM_WORLD,
792  "reduction delta_time = %6.4e delta_time = %6.4e\n",
793  reduction,delta_time
794  ); CHKERRQ(ierr);
795  delta_time *= reduction;
796  if(reduction>1 && delta_time < min_reduction*delta_time0) {
797  delta_time = min_reduction*delta_time0;
798  }
799 
800  ierr = DMoFEMMeshToGlobalVector(
801  dm,D,INSERT_VALUES,SCATTER_REVERSE
802  ); CHKERRQ(ierr);
803  ierr = DMoFEMLoopFiniteElements(
804  dm,"PLASTIC",&small_strain_plasticity.feUpdate
805  ); CHKERRQ(ierr);
806 
807 
808  //Save data on mesh
809  ierr = DMoFEMLoopFiniteElements(
810  dm,"PLASTIC",&post_proc
811  ); CHKERRQ(ierr);
812  string out_file_name;
813  {
814  std::ostringstream stm;
815  stm << "out_" << step << ".h5m";
816  out_file_name = stm.str();
817  }
818  ierr = PetscPrintf(
819  PETSC_COMM_WORLD,"out file %s\n",out_file_name.c_str()
820  ); CHKERRQ(ierr);
821  rval = post_proc.postProcMesh.write_file(
822  out_file_name.c_str(),"MOAB","PARALLEL=WRITE_PART"
823  ); CHKERRQ_MOAB(rval);
824 
825  //===================================== Homgenised stress (for given strain) ====================================================
826  VectorDouble StressHomo;
827  StressHomo.resize(6);
828  StressHomo.clear();
829 
830  //calculate honmogenised stress
831  //create a vector for 6 components of homogenized stress
832  Vec stress_homo;
833  int stress_homo_ghost[] = { 0,1,2,3,4,5,6 };
834  ierr = VecCreateGhost(
835  PETSC_COMM_WORLD,(!m_field.get_comm_rank())?6:0,6,6,stress_homo_ghost,&stress_homo
836  ); CHKERRQ(ierr);
837 
838  lagrangian_element_periodic.setRVEBCsHomoStressOperators("DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",stress_homo);
839 
840  PetscScalar *avec;
841  ierr = VecGetArray(stress_homo,&avec); CHKERRQ(ierr);
842  ierr = VecZeroEntries(stress_homo); CHKERRQ(ierr);
843  ierr = m_field.loop_finite_elements(
844  "PLASTIC_PROB","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCStress()
845  ); CHKERRQ(ierr);
846  ierr = PetscOptionsGetReal(
847  PETSC_NULL,"-my_rve_volume",&rve_volume,PETSC_NULL
848  ); CHKERRQ(ierr);
849  ierr = VecAssemblyBegin(stress_homo); CHKERRQ(ierr);
850  ierr = VecAssemblyEnd(stress_homo); CHKERRQ(ierr);
851  ierr = VecScale(stress_homo,1.0/rve_volume); CHKERRQ(ierr);
852  ierr = VecGhostUpdateBegin(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
853  ierr = VecGhostUpdateEnd(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
854 
855  for(int jj=0; jj<6; jj++) {
856  StressHomo(jj) = avec[jj];
857  }
858 
859  double scale;
860  ierr = time_force_scale.getForceScale(t,scale); CHKERRQ(ierr);
861 
862  PetscPrintf(PETSC_COMM_WORLD,
863  "Macro_Strain Homo_Stress Path "
864  );
865 
866  //cout<<"Macro-strain ";
867  for(int ii=0; ii<6; ii++) {
868  PetscPrintf(
869  PETSC_COMM_WORLD,
870  "%8.5e\t",
871  t*given_strain(ii)
872  );
873  }
874 
875  //cout<<"Homo stress ";
876  for(int ii=0; ii<6; ii++) {
877  PetscPrintf(
878  PETSC_COMM_WORLD,
879  "%8.5e\t",
880  StressHomo(ii)
881  );
882  }
883 
884  PetscPrintf(PETSC_COMM_WORLD,
885  "\n"
886  );
887 
888  //==============================================================================================================================
889 
890  }
891  }
892 
893  ierr = VecDestroy(&D0); CHKERRQ(ierr);
894  ierr = MatDestroy(&Aij); CHKERRQ(ierr);
895  ierr = VecDestroy(&F); CHKERRQ(ierr);
896  ierr = VecDestroy(&D); CHKERRQ(ierr);
897 
898  }
899  CATCH_ERRORS;
900 
901  MoFEM::Core::Finalize();
902 }
SmallStrainPlasticity::OpGetCommonDataAtGaussPts
Definition: SmallStrainPlasticity.hpp:413
MethodForForceScaling.hpp
SIDESET
@ SIDESET
Definition: definitions.h:147
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:215
SmallStrainParaboloidalPlasticity::nt
double nt
Definition: SmallStrainPlasticityMaterialModels.hpp:175
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
BCs_RVELagrange_Periodic
Definition: BCs_RVELagrange_Periodic.hpp:19
SmallStrainPlasticity::feUpdate
MyVolumeFE feUpdate
Definition: SmallStrainPlasticity.hpp:362
SmallStrainParaboloidalPlasticity::lambda
double lambda
Definition: SmallStrainPlasticityMaterialModels.hpp:171
BCs_RVELagrange_Disp.hpp
EntityHandle
Face_CenPos_Handle_multiIndex
multi_index_container< Face_CenPos_Handle, indexed_by< ordered_non_unique< tag< xcoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::xcoord > >, ordered_non_unique< tag< ycoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::ycoord > >, ordered_non_unique< tag< zcoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::zcoord > >, ordered_unique< tag< Tri_Hand_tag >, member< Face_CenPos_Handle, const EntityHandle,&Face_CenPos_Handle::Tri_Hand > >, ordered_unique< tag< Composite_xyzcoord >, composite_key< Face_CenPos_Handle, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::xcoord >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::ycoord >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::zcoord > > >> > Face_CenPos_Handle_multiIndex
Definition: homogenisation_periodic_atom_test.cpp:96
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:172
SmallStrainPlasticity::OpUpdate
Definition: SmallStrainPlasticity.hpp:418
PostProcTemplateOnRefineMesh::postProcMesh
moab::Interface & postProcMesh
Definition: PostProcOnRefMesh.hpp:122
BCs_RVELagrange_Periodic::setRVEBCsOperatorsNonlinear
PetscErrorCode setRVEBCsOperatorsNonlinear(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &fvec, Vec f, VectorDouble given_strain)
Definition: BCs_RVELagrange_Periodic.hpp:1091
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:454
SmallStrainParaboloidalPlasticity
Small strain plasticity with paraboloidal yield criterion (Isotropic Hardening)
Definition: SmallStrainPlasticityMaterialModels.hpp:162
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
SmallStrainPlasticity::ClosestPointProjection::plasticStrain
VectorDouble plasticStrain
Definition: SmallStrainPlasticity.hpp:489
BCs_RVELagrange_Periodic::methodsOp
boost::ptr_vector< MethodForForceScaling > methodsOp
Definition: BCs_RVELagrange_Periodic.hpp:21
SmallStrainPlasticity::ClosestPointProjection::sTrain
VectorDouble sTrain
Definition: SmallStrainPlasticity.hpp:485
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
LAMBDA
#define LAMBDA(E, NU)
Definition: fem_tools.h:22
BCs_RVELagrange_Periodic::feRVEBCLhs
MyPrismFE feRVEBCLhs
Definition: BCs_RVELagrange_Periodic.hpp:34
MoFEM.hpp
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:527
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:541
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
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_mech_plas_periodic.cpp:107
SmallStrainPlasticity::ClosestPointProjection::tAgs
vector< int > tAgs
Definition: SmallStrainPlasticity.hpp:506
out_file_name
char out_file_name[255]
Definition: initial_diffusion.cpp:53
SmallStrainPlasticity::ClosestPointProjection::snesCreate
PetscErrorCode snesCreate()
Definition: SmallStrainPlasticity.cpp:1027
help
static char help[]
Definition: rve_mech_plas_periodic.cpp:61
SmallStrainPlasticity::ClosestPointProjection::createMatAVecR
virtual PetscErrorCode createMatAVecR()
Definition: SmallStrainPlasticity.cpp:924
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::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:501
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
BCs_RVELagrange_Periodic::feRVEBCRhs
MyPrismFE feRVEBCRhs
Definition: BCs_RVELagrange_Periodic.hpp:33
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
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
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1171
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
SmallStrainParaboloidalPlasticity::sIgma_yt
double sIgma_yt
Definition: SmallStrainPlasticityMaterialModels.hpp:174
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1201
xcoord_tag
Definition: homogenisation_periodic_atom_test.cpp:68
ContactOps::scale
double scale
Definition: EshelbianContact.hpp:22
SmallStrainPlasticity::ClosestPointProjection::tOl
double tOl
Definition: SmallStrainPlasticity.hpp:491
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
SmallStrainParaboloidalPlasticity::nc
double nc
Definition: SmallStrainPlasticityMaterialModels.hpp:175
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
TimeForceScale
Force scale operator for reading two columns.
Definition: TimeForceScale.hpp:18
MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:152
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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.
double
BCs_RVELagrange_Trac_Rigid_Trans.hpp
SmallStrainPlasticity::feLhs
MyVolumeFE feLhs
Definition: SmallStrainPlasticity.hpp:361
TimeForceScale.hpp
MoFEM::DMoFEMMeshToGlobalVector
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:539
MoFEM::DMMoFEMGetSnesCtx
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMoFEM.cpp:1098
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
SmallStrainPlasticityMaterialModels.hpp
Small Strain plasticity material models.
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:173
SmallStrainPlasticity::commonData
CommonData commonData
Definition: SmallStrainPlasticity.hpp:394
Projection10NodeCoordsOnField.hpp
BCs_RVELagrange_Trac
Definition: BCs_RVELagrange_Trac.hpp:22
MoFEM::DMMoFEMSNESSetJacobian
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Definition: DMMoFEM.cpp:763
PostProcTemplateVolumeOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Definition: PostProcOnRefMesh.hpp:301
PostProcTemplateOnRefineMesh::addFieldValuesGradientPostProc
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
Definition: PostProcOnRefMesh.hpp:195
PostProcOnRefMesh.hpp
Post-process fields on refined mesh.
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
BCs_RVELagrange_Periodic::getLoopFeRVEBCRhs
MyPrismFE & getLoopFeRVEBCRhs()
Definition: BCs_RVELagrange_Periodic.hpp:39
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
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
SmallStrainPlasticity
Small strain finite element implementation.
Definition: SmallStrainPlasticity.hpp:320
MoFEM::DMMoFEMCreateMoFEM
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:118
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
SmallStrainPlasticity.hpp
Operators and data structures for small strain plasticity.
Face_CenPos_Handle::Face_CenPos_Handle
Face_CenPos_Handle(double _xcoord, double _ycoord, double _zcoord, const EntityHandle _Tri_Hand)
Definition: rve_mech_plas_periodic.cpp:75
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:23
SmallStrainPlasticity::ClosestPointProjection::internalVariables
VectorDouble internalVariables
Definition: SmallStrainPlasticity.hpp:486
SmallStrainParaboloidalPlasticity::nup
double nup
Definition: SmallStrainPlasticityMaterialModels.hpp:172
MoFEM::SnesRhs
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition: SnesCtx.cpp:27
Range
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition: DeprecatedPetsc.hpp:192
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
TimeForceScale::getForceScale
MoFEMErrorCode getForceScale(const double ts_t, double &scale)
Definition: TimeForceScale.hpp:89
SmallStrainParaboloidalPlasticity::mu
double mu
Definition: SmallStrainPlasticityMaterialModels.hpp:170
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:955
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
std
Definition: enable_if.hpp:5
BCs_RVELagrange_Periodic.hpp
SmallStrainParaboloidalPlasticity::Hc
double Hc
Definition: SmallStrainPlasticityMaterialModels.hpp:173
SmallStrainPlasticity::OpAssembleRhs
Definition: SmallStrainPlasticity.hpp:438
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
SmallStrainPlasticity::createInternalVariablesTag
PetscErrorCode createInternalVariablesTag()
Definition: SmallStrainPlasticity.cpp:34
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
SmallStrainParaboloidalPlasticity::Ht
double Ht
Definition: SmallStrainPlasticityMaterialModels.hpp:173
Tri_Hand_tag
Definition: homogenisation_periodic_atom_test.cpp:71
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::DeprecatedCoreInterface::seed_ref_level_3D
DEPRECATED MoFEMErrorCode seed_ref_level_3D(const EntityHandle meshset, const BitRefLevel &bit, int verb=-1)
seed 2D entities in the meshset and their adjacencies (only TETs adjacencies) in a particular BitRefL...
Definition: DeprecatedCoreInterface.cpp:18
MoFEM::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
SmallStrainParaboloidalPlasticity::sIgma_yc
double sIgma_yc
Definition: SmallStrainPlasticityMaterialModels.hpp:174
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
SmallStrainPlasticity::CommonData::bBar
bool bBar
Definition: SmallStrainPlasticity.hpp:389
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
SmallStrainPlasticity::OpAssembleLhs
Definition: SmallStrainPlasticity.hpp:456
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
SmallStrainPlasticity::feRhs
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
Definition: SmallStrainPlasticity.hpp:360
main
int main(int argc, char *argv[])
Definition: rve_mech_plas_periodic.cpp:116
MoFEM::SnesCtx
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:13
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
MoFEM::SnesMat
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
Definition: SnesCtx.cpp:139
BCs_RVELagrange_Trac.hpp
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:590
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
ycoord_tag
Definition: homogenisation_periodic_atom_test.cpp:69
BCs_RVELagrange_Periodic::getLoopFeRVEBCStress
MyPrismFE & getLoopFeRVEBCStress()
Definition: BCs_RVELagrange_Periodic.hpp:42
BCs_RVELagrange_Disp::feRVEBCLhs
MyTriFE feRVEBCLhs
Definition: BCs_RVELagrange_Disp.hpp:30
MU
#define MU(E, NU)
Definition: fem_tools.h:23
F
@ F
Definition: free_surface.cpp:394
BCs_RVELagrange_Periodic::feRVEBCRhsResidual
MyPrismFE feRVEBCRhsResidual
Definition: BCs_RVELagrange_Periodic.hpp:36
NOFIELD
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1127
MoFEM::DMMoFEMSNESSetFunction
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Definition: DMMoFEM.cpp:722
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
PostProcTemplateOnRefineMesh::addFieldValuesPostProc
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
Definition: PostProcOnRefMesh.hpp:153
zcoord_tag
Definition: homogenisation_periodic_atom_test.cpp:70
SmallStrainPlasticity::OpCalculateStress
Definition: SmallStrainPlasticity.hpp:603
Generated by Doxygen 1.8.17 and hosted at University of Glasgow