v0.13.1
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>
23using namespace MoFEM;
24
25#include <boost/program_options.hpp>
26using namespace std;
27namespace 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
36#include <TimeForceScale.hpp>
37#include <VolumeCalculation.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>
50
52#include <PostProcOnRefMesh.hpp>
53
54#include <string>
55
56using namespace boost::numeric;
57
58
59
60
61static 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
72struct 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
79struct xcoord_tag {}; //tags to used in the multindex container
80struct ycoord_tag {};
81struct zcoord_tag {};
82struct Tri_Hand_tag {};
83struct Composite_xyzcoord {};
84
85typedef multi_index_container<
87indexed_by<
88ordered_non_unique<
89tag<xcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord> >,
90
91ordered_non_unique<
92tag<ycoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord> >,
93
94ordered_non_unique<
95tag<zcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> >,
96
97ordered_unique<
98tag<Tri_Hand_tag>, member<Face_CenPos_Handle,const EntityHandle,&Face_CenPos_Handle::Tri_Hand> >,
99
100ordered_unique<
101tag<Composite_xyzcoord>,
102composite_key<
104member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord>,
105member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord>,
106member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> > >
108
109//============================================================================================
110//============================================================================================
111
112
113
114
115
116int 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
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);
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(
616 "DISPLACEMENT",small_strain_plasticity.commonData
617 )
618 );
619 small_strain_plasticity.feRhs.getOpPtrVector().push_back(
621 m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp,false
622 )
623 );
624 small_strain_plasticity.feRhs.getOpPtrVector().push_back(
626 "DISPLACEMENT",small_strain_plasticity.commonData
627 )
628 );
629 small_strain_plasticity.feLhs.getOpPtrVector().push_back(
631 "DISPLACEMENT",small_strain_plasticity.commonData
632 )
633 );
634 small_strain_plasticity.feLhs.getOpPtrVector().push_back(
636 m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp,false
637 )
638 );
639 small_strain_plasticity.feLhs.getOpPtrVector().push_back(
641 "DISPLACEMENT",small_strain_plasticity.commonData
642 )
643 );
644 small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
646 "DISPLACEMENT",small_strain_plasticity.commonData
647 )
648 );
649 small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
651 m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp,false
652 )
653 );
654 small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
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
801 dm,D,INSERT_VALUES,SCATTER_REVERSE
802 ); CHKERRQ(ierr);
804 dm,"PLASTIC",&small_strain_plasticity.feUpdate
805 ); CHKERRQ(ierr);
806
807
808 //Save data on mesh
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);
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 }
900
902}
Post-process fields on refined mesh.
Operators and data structures for small strain plasticity.
Small Strain plasticity material models.
Operator can be used with any volume element to calculate sum of volumes of all volumes in the set.
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ NOBASE
Definition: definitions.h:72
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:554
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:97
@ H1
continuous field
Definition: definitions.h:98
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
@ SIDESET
Definition: definitions.h:160
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:467
#define MU(E, NU)
Definition: fem_tools.h:33
#define LAMBDA(E, NU)
Definition: fem_tools.h:32
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1082
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: DMMMoFEM.cpp:677
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: DMMMoFEM.cpp:130
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:462
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:482
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: DMMMoFEM.cpp:718
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMMoFEM.cpp:1156
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:545
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMMoFEM.cpp:1126
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMMoFEM.cpp:1053
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:494
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
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
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
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
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.
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
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.
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.
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
char mesh_file_name[255]
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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:148
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
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:39
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
CoreTmp< 0 > Core
Definition: Core.hpp:1096
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
const double D
diffusivity
STL namespace.
double young_modulus
Definition: plastic.cpp:105
double scale
Definition: plastic.cpp:103
constexpr double t
plate stiffness
Definition: plate.cpp:76
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
int main(int argc, char *argv[])
static char help[]
map< int, RVEBC_Data > setOfRVEBC
maps side set id with appropriate FluxData
PetscErrorCode addLagrangiangElement(const string element_name, const string field_name, const string lagrang_field_name, const string mesh_nodals_positions)
boost::ptr_vector< MethodForForceScaling > methodsOp
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec stress_homo)
PetscErrorCode addLagrangiangElement(const string element_name, const string field_name, const string lagrang_field_name, const string mesh_nodals_positions, Range &periodic_prisms)
PetscErrorCode setRVEBCsOperatorsNonlinear(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &fvec, Vec f, VectorDouble given_strain)
PetscErrorCode setRVEBCsRigidBodyTranOperators(string field_name, string lagrang_field_name, Mat aij, map< int, RVEBC_Data > setOfRVEBC)
Face_CenPos_Handle(double _xcoord, double _ycoord, double _zcoord, const EntityHandle _Tri_Hand)
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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.
virtual int get_comm_rank() const =0
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
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...
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Projection of edge entities with one mid-node on hierarchical basis.
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:27
moab::Interface & postProcMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Post processing.
Small strain plasticity with paraboloidal yield criterion (Isotropic Hardening)
Small strain finite element implementation.
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
PetscErrorCode createInternalVariablesTag()
Force scale operator for reading two columns.
MoFEMErrorCode getForceScale(const double ts_t, double &scale)
Calculate volume.