v0.14.0
Loading...
Searching...
No Matches
Classes | Typedefs | Functions | Variables
rve_mech_plas_periodic.cpp File Reference
#include <MoFEM.hpp>
#include <boost/program_options.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <MethodForForceScaling.hpp>
#include <TimeForceScale.hpp>
#include <VolumeCalculation.hpp>
#include <BCs_RVELagrange_Disp.hpp>
#include <BCs_RVELagrange_Trac.hpp>
#include "BCs_RVELagrange_Periodic.hpp"
#include <BCs_RVELagrange_Trac_Rigid_Trans.hpp>
#include <adolc/adolc.h>
#include <SmallStrainPlasticity.hpp>
#include <SmallStrainPlasticityMaterialModels.hpp>
#include <Projection10NodeCoordsOnField.hpp>
#include <PostProcOnRefMesh.hpp>
#include <string>

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
 

Functions

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

Variables

static char help []
 

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 107 of file rve_mech_plas_periodic.cpp.

Function Documentation

◆ main()

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

Definition at line 116 of file rve_mech_plas_periodic.cpp.

116 {
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}
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition definitions.h:385
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
NOBASE
@ NOBASE
Definition definitions.h:59
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:84
H1
@ H1
continuous field
Definition definitions.h:85
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition definitions.h:228
SIDESET
@ SIDESET
Definition definitions.h:160
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition definitions.h:467
order
constexpr int order
Definition dg_projection.cpp:18
MU
#define MU(E, NU)
Definition fem_tools.h:23
LAMBDA
#define LAMBDA(E, NU)
Definition fem_tools.h:22
F
@ F
Definition free_surface.cpp:394
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition DMMoFEM.cpp:1123
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:497
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:718
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:114
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:523
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:759
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition DMMoFEM.cpp:1197
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
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:586
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition DMMoFEM.cpp:1167
MoFEM::DMMoFEMGetSnesCtx
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition DMMoFEM.cpp:1094
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:535
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::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::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_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::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
MoFEM::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::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::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::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.
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
out_file_name
char out_file_name[255]
Definition initial_diffusion.cpp:53
D
double D
Definition initial_diffusion.cpp:44
mesh_file_name
char mesh_file_name[255]
Definition mesh_smoothing.cpp:22
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition Exceptions.hpp:76
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition Exceptions.hpp:74
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
MoFEM::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
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition DeprecatedPetsc.hpp:142
MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition DeprecatedPetsc.hpp:152
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
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition DeprecatedPetsc.hpp:182
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition DeprecatedPetsc.hpp:192
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition DeprecatedPetsc.hpp:172
young_modulus
double young_modulus
Young modulus.
Definition plastic.cpp:121
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition plastic.cpp:122
scale
double scale
Definition plastic.cpp:119
t
constexpr double t
plate stiffness
Definition plate.cpp:58
help
static char help[]
Definition rve_mech_plas_periodic.cpp:61
BCs_RVELagrange_Periodic
Definition BCs_RVELagrange_Periodic.hpp:19
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
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_comm_rank
virtual int get_comm_rank() const =0
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition Core.hpp:82
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
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
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::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition Projection10NodeCoordsOnField.hpp:24
MoFEM::SnesCtx
Interface for nonlinear (SNES) solver.
Definition SnesCtx.hpp:13
PostProcVolumeOnRefinedMesh
Post processing.
Definition PostProcOnRefMesh.hpp:957
TimeForceScale
Force scale operator for reading two columns.
Definition TimeForceScale.hpp:18
Tri_Hand_tag
Definition homogenisation_periodic_atom_test.cpp:71
VolumeCalculation
Calculate volume.
Definition VolumeCalculation.hpp:18
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

Variable Documentation

◆ help

char help[]
static
Initial value:
=
"...\n"
"\n"

Definition at line 61 of file rve_mech_plas_periodic.cpp.

Generated by Doxygen 1.12.0 and hosted at University of Glasgow