v0.13.1
rve_mech_plas_interface_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
35
36#include <MethodForForceScaling.hpp>
37#include <TimeForceScale.hpp>
38#include <VolumeCalculation.hpp>
39
40#include <VolumeCalculation.hpp>
41#include <BCs_RVELagrange_Disp.hpp>
42#include <BCs_RVELagrange_Trac.hpp>
43#include "BCs_RVELagrange_Periodic.hpp"
44#include <BCs_RVELagrange_Trac_Rigid_Trans.hpp>
45
46#include <CohesiveInterfaceElement.hpp>
47using namespace CohesiveElement;
48
49#ifndef WITH_ADOL_C
50#error "MoFEM need to be compiled with ADOL-C"
51#endif
52#include <adolc/adolc.h>
53#include <NonLinearElasticElement.hpp>
54#include <SmallTransverselyIsotropic.hpp>
55
56#include <adolc/adolc.h>
57#include <SmallStrainPlasticity.hpp>
58#include <SmallStrainPlasticityMaterialModels.hpp>
59
60
61#include <Projection10NodeCoordsOnField.hpp>
62#include <PostProcOnRefMesh.hpp>
63
64#include <string>
65
66using namespace boost::numeric;
67
68
69
70
71static char help[] =
72"...\n"
73"\n";
74
75
76//=================================================================================================================================
77//Define class and multindex container to store data for traiangles on the boundary of the RVE (it cannot be defined within main)
78//=================================================================================================================================
79
80struct Face_CenPos_Handle {
81 double xcoord, ycoord, zcoord;
82 const EntityHandle Tri_Hand;
83 Face_CenPos_Handle(double _xcoord, double _ycoord, double _zcoord, const EntityHandle _Tri_Hand):xcoord(_xcoord),
84 ycoord(_ycoord), zcoord(_zcoord), Tri_Hand(_Tri_Hand) {}
85};
86
87struct xcoord_tag {}; //tags to used in the multindex container
88struct ycoord_tag {};
89struct zcoord_tag {};
90struct Tri_Hand_tag {};
91struct Composite_xyzcoord {};
92
93typedef multi_index_container<
94Face_CenPos_Handle,
95indexed_by<
96ordered_non_unique<
97tag<xcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord> >,
98
99ordered_non_unique<
100tag<ycoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord> >,
101
102ordered_non_unique<
103tag<zcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> >,
104
105ordered_unique<
106tag<Tri_Hand_tag>, member<Face_CenPos_Handle,const EntityHandle,&Face_CenPos_Handle::Tri_Hand> >,
107
108ordered_unique<
109tag<Composite_xyzcoord>,
110composite_key<
111Face_CenPos_Handle,
112member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord>,
113member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord>,
114member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> > >
115> > Face_CenPos_Handle_multiIndex;
116
117//============================================================================================
118//============================================================================================
119
120
121
122int main(int argc, char *argv[]) {
123 MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
124
125 try {
126
127 moab::Core mb_instance;
128 moab::Interface& moab = mb_instance;
129 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
130 if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
131
132 PetscBool flg = PETSC_TRUE;
133 char mesh_file_name[255];
134 ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
135 if(flg != PETSC_TRUE) {
136 SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
137 }
138
139 PetscInt order;
140 ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
141 if(flg != PETSC_TRUE) {
142 order = 2;
143 }
144 PetscInt bubble_order;
145 ierr = PetscOptionsGetInt(PETSC_NULL,"-my_bubble_order",&bubble_order,&flg); CHKERRQ(ierr);
146 if(flg != PETSC_TRUE) {
147 bubble_order = order;
148 }
149
150 // use this if your mesh is partitioned and you run code on parts,
151 // you can solve very big problems
152 PetscBool is_partitioned = PETSC_FALSE;
153 ierr = PetscOptionsGetBool(PETSC_NULL,"-my_is_partitioned",&is_partitioned,&flg); CHKERRQ(ierr);
154
155 //Applied strain on the RVE (vector of length 6) strain=[xx, yy, zz, xy, xz, zy]^T
156 double mygiven_strain[6];
157 int nmax=6;
158 ierr = PetscOptionsGetRealArray(PETSC_NULL,"-my_given_strain",mygiven_strain,&nmax,&flg); CHKERRQ(ierr);
159 VectorDouble given_strain;
160 given_strain.resize(6);
161 cblas_dcopy(6, &mygiven_strain[0], 1, &given_strain(0), 1);
162 cout<<"given_strain ="<<given_strain<<endl;
163
164 //Read mesh to MOAB
165
166 if(is_partitioned == PETSC_TRUE) {
167 const char *option;
168 option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION=PARALLEL_PARTITION;";
169 rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
170 } else {
171 const char *option;
172 option = "";
173 rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
174 }
175
176
177
178 double young_modulus = 1;
179 SmallStrainParaboloidalPlasticity cp;
180 {
181 cp.tAgs.resize(3);
182 cp.tAgs[0] = 3;
183 cp.tAgs[1] = 4;
184 cp.tAgs[2] = 5;
185 cp.tOl = 1e-12;
186
187 double poisson_ratio = 0.25;
188 cp.sIgma_yt = 1;
189 cp.sIgma_yc = 1;
190
191 cp.Ht = 0.1;
192 cp.Hc = 0.1;
193
194 cp.nup = 0.3;
195
196
197 {
198 ierr = PetscOptionsGetReal(0,"-my_young_modulus",&young_modulus,0); CHKERRQ(ierr);
199 ierr = PetscOptionsGetReal(0,"-my_poisson_ratio",&poisson_ratio,0); CHKERRQ(ierr);
200 cp.mu = MU(young_modulus,poisson_ratio);
201 cp.lambda = LAMBDA(young_modulus,poisson_ratio);
202 ierr = PetscOptionsGetReal(0,"-my_sigma_yt",&cp.sIgma_yt,0); CHKERRQ(ierr);
203 ierr = PetscOptionsGetReal(0,"-my_sigma_yc",&cp.sIgma_yc,0); CHKERRQ(ierr);
204
205 ierr = PetscOptionsGetReal(0,"-my_Ht",&cp.Ht,0); CHKERRQ(ierr);
206 ierr = PetscOptionsGetReal(0,"-my_Hc",&cp.Hc,0); CHKERRQ(ierr);
207
208 ierr = PetscOptionsGetReal(0,"-my_nt",&cp.nt,0); CHKERRQ(ierr);
209 ierr = PetscOptionsGetReal(0,"-my_nc",&cp.nc,0); CHKERRQ(ierr);
210
211
212 ierr = PetscOptionsGetReal(0,"-my_nup",&cp.nup,0); CHKERRQ(ierr);
213 }
214
215 PetscPrintf(PETSC_COMM_WORLD,"young_modulus = %4.2e \n",young_modulus);
216 PetscPrintf(PETSC_COMM_WORLD,"poisson_ratio = %4.2e \n",poisson_ratio);
217
218 PetscPrintf(PETSC_COMM_WORLD,"sIgma_yt = %4.2e \n",cp.sIgma_yt);
219 PetscPrintf(PETSC_COMM_WORLD,"sIgma_yc = %4.2e \n",cp.sIgma_yc);
220
221 PetscPrintf(PETSC_COMM_WORLD,"nup = %4.2e \n",cp.nup);
222
223 cp.sTrain.resize(6,false);
224 cp.sTrain.clear();
225 cp.plasticStrain.resize(6,false);
226 cp.plasticStrain.clear();
227 cp.internalVariables.resize(2,false);
228 cp.internalVariables.clear();
229 cp.createMatAVecR();
230 cp.snesCreate();
231 // ierr = SNESSetFromOptions(cp.sNes); CHKERRQ(ierr);
232 }
233
234
235 // double young_modulus = 1;
236 // SmallStrainJ2Plasticity cp;
237 // {
238 // cp.tAgs.resize(3);
239 // cp.tAgs[0] = 3;
240 // cp.tAgs[1] = 4;
241 // cp.tAgs[2] = 5;
242 // cp.tOl = 1e-12;
243 //
244 // double poisson_ratio = 0.25;
245 // cp.sIgma_y = 1;
246 // cp.H = 0.1;
247 // cp.K = 0;
248 // cp.phi = 1;
249 // {
250 // ierr = PetscOptionsGetReal(0,"-my_young_modulus",&young_modulus,0); CHKERRQ(ierr);
251 // ierr = PetscOptionsGetReal(0,"-my_poisson_ratio",&poisson_ratio,0); CHKERRQ(ierr);
252 // cp.mu = MU(young_modulus,poisson_ratio);
253 // cp.lambda = LAMBDA(young_modulus,poisson_ratio);
254 // ierr = PetscOptionsGetReal(0,"-my_sigma_y",&cp.sIgma_y,0); CHKERRQ(ierr);
255 // ierr = PetscOptionsGetReal(0,"-my_H",&cp.H,0); CHKERRQ(ierr);
256 // ierr = PetscOptionsGetReal(0,"-my_K",&cp.K,0); CHKERRQ(ierr);
257 // ierr = PetscOptionsGetReal(0,"-my_phi",&cp.phi,0); CHKERRQ(ierr);
258 // }
259 //
260 // PetscPrintf(PETSC_COMM_WORLD,"young_modulus = %4.2e \n",young_modulus);
261 // PetscPrintf(PETSC_COMM_WORLD,"poisson_ratio = %4.2e \n",poisson_ratio);
262 // PetscPrintf(PETSC_COMM_WORLD,"sIgma_y = %4.2e \n",cp.sIgma_y);
263 // PetscPrintf(PETSC_COMM_WORLD,"H = %4.2e \n",cp.H);
264 // PetscPrintf(PETSC_COMM_WORLD,"K = %4.2e \n",cp.K);
265 // PetscPrintf(PETSC_COMM_WORLD,"phi = %4.2e \n",cp.phi);
266 //
267 //
268 // cp.sTrain.resize(6,false);
269 // cp.sTrain.clear();
270 // cp.plasticStrain.resize(6,false);
271 // cp.plasticStrain.clear();
272 // cp.internalVariables.resize(7,false);
273 // cp.internalVariables.clear();
274 // cp.createMatAVecR();
275 // cp.snesCreate();
276 // // ierr = SNESSetFromOptions(cp.sNes); CHKERRQ(ierr);
277 // }
278
279
280 //Create MoFEM (Joseph) database
281 MoFEM::Core core(moab);
282 MoFEM::Interface& m_field = core;
283
284 vector<BitRefLevel> bit_levels;
285 Tag th_meshset_info;
286 int def_meshset_info[2] = {0,0};
287 rval = moab.tag_get_handle("MESHSET_INFO",2,MB_TYPE_INTEGER,th_meshset_info,MB_TAG_CREAT|MB_TAG_SPARSE,&def_meshset_info);
288 int meshset_data[2];
289 EntityHandle root = moab.get_root_set();
290 rval = moab.tag_get_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
291 if(meshset_data[0]==0) {
292 meshset_data[0] = 1;
293 rval = moab.tag_set_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
294 }
295 bit_levels.push_back(BitRefLevel().set(meshset_data[0]-1));
296
297
298// EntityHandle out_meshset;
299// rval = moab.create_meshset(MESHSET_SET,out_meshset); CHKERRQ_MOAB(rval);
300// // ierr = m_field.get_problem_finite_elements_entities("POTENTIAL_PROBLEM","POTENTIAL_ELEM",out_meshset); CHKERRQ(ierr);
301// ierr = m_field.get_entities_by_ref_level(bit_levels.back(),BitRefLevel().set(),out_meshset); CHKERRQ(ierr);
302// Range LatestRefinedTets;
303// rval = moab.get_entities_by_type(out_meshset, MBTET,LatestRefinedTets,true); CHKERRQ_MOAB(rval);
304// Range LatestRefinedPrisms;
305// rval = moab.get_entities_by_type(out_meshset, MBPRISM,LatestRefinedPrisms,true); CHKERRQ_MOAB(rval);
306
307
308 BitRefLevel problem_bit_level = bit_levels.back();
309
310 // const clock_t begin_time = clock();
311 ierr = m_field.build_fields(); CHKERRQ(ierr);
312 ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
313
314
315 Range interface_prims_on_problem_bit_level;
316 ierr = m_field.get_entities_by_type_and_ref_level(
317 problem_bit_level,BitRefLevel().set(),MBPRISM,interface_prims_on_problem_bit_level
318 ); CHKERRQ(ierr);
319 Range tets_on_problem_bit_level;
320 ierr = m_field.get_entities_by_type_and_ref_level(
321 problem_bit_level,BitRefLevel().set(),MBTET,tets_on_problem_bit_level
322 ); CHKERRQ(ierr);
323
324 cout<<"interface_prims_on_problem_bit_level "<<interface_prims_on_problem_bit_level.size()<<endl;
325 cout<<"tets_on_problem_bit_level "<<tets_on_problem_bit_level.size()<<endl;
326
327 //=======================================================================================================
328 //Add Periodic Prisims Between Triangles on -ve and +ve faces to implement periodic bounary conditions
329 //=======================================================================================================
330 Range preriodic_prisms;
331 //Populating the Multi-index container with -ve triangles
332
333
334 Range SurTrisNeg; //Old + new traiangles on side set 101
335 ierr = m_field.get_cubit_msId_entities_by_dimension(101,SIDESET,2,SurTrisNeg,true); CHKERRQ(ierr);
336 ierr = PetscPrintf(PETSC_COMM_WORLD,"Old+New triangles on SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
337
338 Range AllTris; //All new traiangles in the new-bitlevel
339 ierr = m_field.get_entities_by_type_and_ref_level(problem_bit_level,BitRefLevel().set(),MBTRI,AllTris); CHKERRQ(ierr);
340 ierr = PetscPrintf(PETSC_COMM_WORLD,"All triangles in the new bit-level = %d\n",AllTris.size()); CHKERRQ(ierr);
341
342 //Only triangles in the new bit-level
343 SurTrisNeg = intersect(AllTris,SurTrisNeg);
344 ierr = PetscPrintf(PETSC_COMM_WORLD,"New triangles on SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
345
346// cout<<"SurTrisNeg "<<SurTrisNeg<<endl;
347 Face_CenPos_Handle_multiIndex Face_CenPos_Handle_varNeg, Face_CenPos_Handle_varPos;
348 double TriCen[3], coords_Tri[9];
349
350
351
352 double roundfact=1e3;
353
354 for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
355 const EntityHandle* conn_face; int num_nodes_Tri;
356
357 //get nodes attached to the face
358 rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
359 //get nodal coordinates
360 rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
361
362 //Find triangle centriod
363 TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
364 TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
365 TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
366
367// cout<<"\nTriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
368 TriCen[0]=round(TriCen[0]*roundfact)/roundfact;
369 TriCen[1]=round(TriCen[1]*roundfact)/roundfact;
370 TriCen[2]=round(TriCen[2]*roundfact)/roundfact;
371
372// cout<<"TriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
373 //fill the multi-index container with centriod coordinates and triangle handles
374 Face_CenPos_Handle_varNeg.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
375 // for(int ii=0; ii<3; ii++) cout<<"TriCen "<<TriCen[ii]<<endl;
376 // cout<<endl<<endl;
377 }
378
379 // double aaa;
380 // aaa=0.5011;
381 // cout<<"\n\n\n\n\nfloor(aaa+0.5) = "<<double(int(aaa*roundfact+0.5))/roundfact<<endl<<endl<<endl<<endl;
382 // aaa=-0.5011;
383 // cout<<"\n\n\n\n\nfloor(aaa+0.5) = "<<double(int(aaa*roundfact-0.5))/roundfact<<endl<<endl<<endl<<endl;
384
385
386 //Populating the Multi-index container with +ve triangles
387 Range SurTrisPos;
388 ierr = m_field.get_cubit_msId_entities_by_dimension(102,SIDESET,2,SurTrisPos,true); CHKERRQ(ierr);
389 ierr = PetscPrintf(PETSC_COMM_WORLD,"Old+New triangles on SideSet 102 = %d\n",SurTrisPos.size()); CHKERRQ(ierr);
390
391 //Only triangles in the new bit-level
392 SurTrisPos = intersect(AllTris,SurTrisPos);
393 ierr = PetscPrintf(PETSC_COMM_WORLD,"New triangles on SideSet 102 = %d\n",SurTrisPos.size()); CHKERRQ(ierr);
394
395
396
397
398 for(Range::iterator it = SurTrisPos.begin(); it!=SurTrisPos.end(); it++) {
399 const EntityHandle* conn_face; int num_nodes_Tri;
400
401 //get nodes attached to the face
402 rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
403 //get nodal coordinates
404 rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
405
406 //Find triangle centriod
407 TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
408 TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
409 TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
410
411 //round values to 3 disimal places
412 TriCen[0]=round(TriCen[0]*roundfact)/roundfact;
413 TriCen[1]=round(TriCen[1]*roundfact)/roundfact;
414 TriCen[2]=round(TriCen[2]*roundfact)/roundfact;
415
416// cout<<"TriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
417
418 //fill the multi-index container with centriod coordinates and triangle handles
419 Face_CenPos_Handle_varPos.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
420 }
421
422 //Find minimum and maximum X, Y and Z coordinates of the RVE (using multi-index container)
423 double XcoordMin, YcoordMin, ZcoordMin, XcoordMax, YcoordMax, ZcoordMax;
424 typedef Face_CenPos_Handle_multiIndex::index<xcoord_tag>::type::iterator Tri_Xcoord_iterator;
425 typedef Face_CenPos_Handle_multiIndex::index<ycoord_tag>::type::iterator Tri_Ycoord_iterator;
426 typedef Face_CenPos_Handle_multiIndex::index<zcoord_tag>::type::iterator Tri_Zcoord_iterator;
427 Tri_Xcoord_iterator XcoordMin_it, XcoordMax_it;
428 Tri_Ycoord_iterator YcoordMin_it, YcoordMax_it;
429 Tri_Zcoord_iterator ZcoordMin_it, ZcoordMax_it;
430
431 //XcoordMax_it-- because .end() will point iterator after the data range but .begin() will point the iteratore to the first value of range
432 XcoordMin_it=Face_CenPos_Handle_varNeg.get<xcoord_tag>().begin(); XcoordMin=XcoordMin_it->xcoord;
433 XcoordMax_it=Face_CenPos_Handle_varPos.get<xcoord_tag>().end(); XcoordMax_it--; XcoordMax=XcoordMax_it->xcoord;
434 YcoordMin_it=Face_CenPos_Handle_varNeg.get<ycoord_tag>().begin(); YcoordMin=YcoordMin_it->ycoord;
435 YcoordMax_it=Face_CenPos_Handle_varPos.get<ycoord_tag>().end(); YcoordMax_it--; YcoordMax=YcoordMax_it->ycoord;
436 ZcoordMin_it=Face_CenPos_Handle_varNeg.get<zcoord_tag>().begin(); ZcoordMin=ZcoordMin_it->zcoord;
437 ZcoordMax_it=Face_CenPos_Handle_varPos.get<zcoord_tag>().end(); ZcoordMax_it--; ZcoordMax=ZcoordMax_it->zcoord;
438
439 cout<<"XcoordMin "<<XcoordMin << " XcoordMax "<<XcoordMax <<endl;
440 cout<<"YcoordMin "<<YcoordMin << " YcoordMax "<<YcoordMax <<endl;
441 cout<<"ZcoordMin "<<ZcoordMin << " ZcoordMax "<<ZcoordMax <<endl;
442
443 /*double LxRVE, LyRVE, LzRVE;
444 LxRVE=XcoordMax-XcoordMin;
445 LyRVE=YcoordMax-YcoordMin;
446 LzRVE=ZcoordMax-ZcoordMin;*/
447
448 //Creating Prisims between triangles on -ve and +ve faces
449 typedef Face_CenPos_Handle_multiIndex::index<Tri_Hand_tag>::type::iterator Tri_Hand_iterator;
450 Tri_Hand_iterator Tri_Neg;
451 typedef Face_CenPos_Handle_multiIndex::index<Composite_xyzcoord>::type::iterator xyzcoord_iterator;
452 xyzcoord_iterator Tri_Pos;
453 double XPos, YPos, ZPos;
454 //int count=0;
455
456 //loop over -ve triangles to create prisims elemenet between +ve and -ve triangles
457 int count=1;
458 int count1=0;
459 for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
460
461// count1++;
462// cout<<"count1 = "<<count1<< endl;
463// if(count1 == 24) {
464
465
466 Tri_Neg=Face_CenPos_Handle_varNeg.get<Tri_Hand_tag>().find(*it);
467 //cout<<"xcoord= "<<Tri_iit->xcoord << " ycoord= "<< Tri_iit->ycoord << " ycoord= "<< Tri_iit->zcoord <<endl;
468
469 //corresponding +ve triangle
470 if(Tri_Neg->xcoord==XcoordMin){XPos=XcoordMax; YPos=Tri_Neg->ycoord; ZPos=Tri_Neg->zcoord;};
471 if(Tri_Neg->ycoord==YcoordMin){XPos=YPos=Tri_Neg->xcoord; YPos=YcoordMax; ZPos=Tri_Neg->zcoord;};
472 if(Tri_Neg->zcoord==ZcoordMin){XPos=YPos=Tri_Neg->xcoord; YPos=Tri_Neg->ycoord; ZPos=ZcoordMax; };
473 Tri_Pos=Face_CenPos_Handle_varPos.get<Composite_xyzcoord>().find(boost::make_tuple(XPos, YPos, ZPos));
474
475// cout<<"Tri_Neg->xcoord= "<<Tri_Neg->xcoord << " Tri_Neg->ycoord "<< Tri_Neg->ycoord << " Tri_Neg->zcoord= "<< Tri_Neg->zcoord <<endl;
476// cout<<"Tri_Pos->xcoord= "<<Tri_Pos->xcoord << " Tri_Pos->ycoord "<< Tri_Pos->ycoord << " Tri_Pos->zcoord= "<< Tri_Pos->zcoord <<endl;
477
478
479 //+ve and -ve nodes and their coords (+ve and -ve tiangles nodes can have matching problems, which can produce twisted prism)
480 EntityHandle PrismNodes[6];
481 vector<EntityHandle> TriNodesNeg, TriNodesPos;
482 double CoordNodeNeg[9], CoordNodePos[9];
483 rval = moab.get_connectivity(&(Tri_Neg->Tri_Hand),1,TriNodesNeg,true); CHKERRQ_MOAB(rval);
484 rval = moab.get_connectivity(&(Tri_Pos->Tri_Hand),1,TriNodesPos,true); CHKERRQ_MOAB(rval);
485 rval = moab.get_coords(&TriNodesNeg[0],3,CoordNodeNeg); MOAB_THROW(rval);
486 rval = moab.get_coords(&TriNodesPos[0],3,CoordNodePos); MOAB_THROW(rval);
487 for(int ii=0; ii<3; ii++){
488 PrismNodes[ii]=TriNodesNeg[ii];
489 }
490
491// for(int ii=0; ii<3; ii++){
492// cout<<"xcoord= "<<CoordNodeNeg[3*ii] << " ycoord= "<< CoordNodeNeg[3*ii+1] << " zcoord= "<< CoordNodeNeg[3*ii+2] <<endl;
493// }
494// for(int ii=0; ii<3; ii++){
495// cout<<"xcoord= "<<CoordNodePos[3*ii] << " ycoord= "<< CoordNodePos[3*ii+1] << " zcoord= "<< CoordNodePos[3*ii+2] <<endl;
496// }
497
498 //Match exact nodes to each other to avoide the problem of twisted prisms
499 double XNodeNeg, YNodeNeg, ZNodeNeg, XNodePos, YNodePos, ZNodePos;
500 for(int ii=0; ii<3; ii++){
501 if(Tri_Neg->xcoord==XcoordMin){XNodeNeg=XcoordMax; YNodeNeg=CoordNodeNeg[3*ii+1]; ZNodeNeg=CoordNodeNeg[3*ii+2];};
502 if(Tri_Neg->ycoord==YcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=YcoordMax; ZNodeNeg=CoordNodeNeg[3*ii+2];};
503 if(Tri_Neg->zcoord==ZcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=CoordNodeNeg[3*ii+1]; ZNodeNeg=ZcoordMax;};
504 for(int jj=0; jj<3; jj++){
505 //Round nodal coordinates to 3 dicimal places only for comparison
506 XNodeNeg=round(XNodeNeg*roundfact)/roundfact;
507 YNodeNeg=round(YNodeNeg*roundfact)/roundfact;
508 ZNodeNeg=round(ZNodeNeg*roundfact)/roundfact;
509
510 XNodePos=CoordNodePos[3*jj]; YNodePos=CoordNodePos[3*jj+1]; ZNodePos=CoordNodePos[3*jj+2];
511 XNodePos=round(XNodePos*roundfact)/roundfact;
512 YNodePos=round(YNodePos*roundfact)/roundfact;
513 ZNodePos=round(ZNodePos*roundfact)/roundfact;
514
515 if(XNodeNeg==XNodePos && YNodeNeg==YNodePos && ZNodeNeg==ZNodePos){
516 PrismNodes[3+ii]=TriNodesPos[jj];
517 break;
518 }
519 }
520 }
521 //prism nodes and their coordinates
522 double CoordNodesPrisms[18];
523 rval = moab.get_coords(PrismNodes,6,CoordNodesPrisms); MOAB_THROW(rval);
524// for(int ii=0; ii<6; ii++){
525// cout<<"xcoord= "<<CoordNodesPrisms[3*ii] << " ycoord= "<< CoordNodesPrisms[3*ii+1] << " zcoord= "<< CoordNodesPrisms[3*ii+2] <<endl;
526// }
527// cout<<endl<<endl;
528
529 //insertion of individula prism element and its addition to range preriodic_prisms
530 EntityHandle PeriodicPrism;
531 rval = moab.create_element(MBPRISM,PrismNodes,6,PeriodicPrism); CHKERRQ_MOAB(rval);
532 preriodic_prisms.insert(PeriodicPrism);
533
534// //to see individual prisms
535// Range Prism1;
536// Prism1.insert(PeriodicPrism);
537// EntityHandle out_meshset1;
538// rval = moab.create_meshset(MESHSET_SET,out_meshset1); CHKERRQ_MOAB(rval);
539// rval = moab.add_entities(out_meshset1,Prism1); CHKERRQ_MOAB(rval);
540// ostringstream sss;
541// sss << "Prism" << count << ".vtk"; count++;
542// rval = moab.write_file(sss.str().c_str(),"VTK","",&out_meshset1,1); CHKERRQ_MOAB(rval);
543//
544// string aaa;
545// cin>>aaa;
546//
547// }
548
549 }
550
551
552
553
554 //Adding the newely added prisms to the bit-level (which is known as seeding)
555 /*{
556 Range ents;
557 moab.get_adjacencies(preriodic_prisms,1,true,ents,moab::Interface::UNION);
558 moab.get_adjacencies(preriodic_prisms,2,true,ents,moab::Interface::UNION);
559 ierr = m_field.seed_ref_level(ents.subset_by_type(MBEDGE),problem_bit_level); CHKERRQ(ierr);
560 ierr = m_field.seed_ref_level(ents.subset_by_type(MBQUAD),problem_bit_level); CHKERRQ(ierr);
561 }*/
562 ierr = m_field.seed_ref_level(preriodic_prisms,problem_bit_level); CHKERRQ(ierr);
563 ierr = m_field.seed_finite_elements(preriodic_prisms); CHKERRQ(ierr);
564
565
566 //=======================================================================================================
567 //=======================================================================================================
568 //=======================================================================================================
569
570 EntityHandle out_meshset;
571 rval = moab.create_meshset(MESHSET_SET,out_meshset); CHKERRQ_MOAB(rval);
572 ierr = m_field.get_entities_by_ref_level(bit_levels.back(),BitRefLevel().set(),out_meshset); CHKERRQ(ierr);
573
574 Range LatestRefinedTets;
575 rval = moab.get_entities_by_type(out_meshset, MBTET,LatestRefinedTets,true); CHKERRQ_MOAB(rval);
576
577 Range Interface_Periodic_Prisms;
578 rval = moab.get_entities_by_type(out_meshset, MBPRISM,Interface_Periodic_Prisms,true); CHKERRQ_MOAB(rval);
579 cout<<"Interface_Periodic_Prisms " <<Interface_Periodic_Prisms.size()<<endl;
580
581
582 //Fields
583 ierr = m_field.add_field("DISPLACEMENT",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
584 ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
585 ierr = m_field.add_field("LAGRANGE_MUL_PERIODIC",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr); //For lagrange multipliers to control the periodic motion
586 ierr = m_field.add_field("LAGRANGE_MUL_RIGID_TRANS",NOFIELD,NOBASE,3); CHKERRQ(ierr); //To control the rigid body motion (3 Traslations and 3
587
588 //add entitities (by tets) to the field
589 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"DISPLACEMENT"); CHKERRQ(ierr);
590 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
591
592// ierr = PetscPrintf(PETSC_COMM_WORLD,"number of triangle in SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
593 ierr = m_field.add_ents_to_field_by_type(SurTrisNeg,MBTRI,"LAGRANGE_MUL_PERIODIC"); CHKERRQ(ierr);
594
595 //set app. order
596 ierr = m_field.set_field_order(0,MBTET,"DISPLACEMENT",bubble_order); CHKERRQ(ierr);
597 ierr = m_field.set_field_order(0,MBTRI,"DISPLACEMENT",order); CHKERRQ(ierr);
598 ierr = m_field.set_field_order(0,MBEDGE,"DISPLACEMENT",order); CHKERRQ(ierr);
599 ierr = m_field.set_field_order(0,MBVERTEX,"DISPLACEMENT",1); CHKERRQ(ierr);
600
601 ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
602 ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
603 ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
604 ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
605
606 ierr = m_field.set_field_order(0,MBTRI,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
607 ierr = m_field.set_field_order(0,MBEDGE,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
608 ierr = m_field.set_field_order(0,MBVERTEX,"LAGRANGE_MUL_PERIODIC",1); CHKERRQ(ierr);
609
610 //build field
611 ierr = m_field.build_fields(); CHKERRQ(ierr);
612
613 //10 node tets
614 Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
615 ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material,0); CHKERRQ(ierr);
616
617 //FE
618 ierr = m_field.add_finite_element("PLASTIC"); CHKERRQ(ierr);
619 //Define rows/cols and element data
620 ierr = m_field.modify_finite_element_add_field_row("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
621 ierr = m_field.modify_finite_element_add_field_col("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
622 ierr = m_field.modify_finite_element_add_field_data("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
623 ierr = m_field.modify_finite_element_add_field_data("PLASTIC","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
624
625 Range newtets;
626 for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,BLOCKSET|UNKNOWNNAME,it)) {
627 if(it->getName() != "MAT_PLASTIC") continue;
628 EntityHandle meshset = it->getMeshset();
629 int id = it->getMeshsetId();
630
631 rval = m_field.get_moab().get_entities_by_type(
632 meshset,MBTET,newtets,true
633 ); CHKERRQ_MOAB(rval);
634
635 cout<< "========================== newtets "<<newtets.size()<<endl;
636
637 //intersection of new and old tets for plastic
638 newtets = intersect(newtets,LatestRefinedTets);
639 ierr = m_field.seed_finite_elements(newtets); CHKERRQ(ierr);
640 }
641
642 ierr = m_field.add_ents_to_finite_element_by_type(newtets,MBTET,"PLASTIC"); CHKERRQ(ierr);
643 cout<< "========================== newtets "<<newtets.size()<<endl;
644
645
646 //================================================================================================================================
647 // Trans-Isotropic Yarns
648 //================================================================================================================================
649
650 NonlinearElasticElement trans_elastic(m_field,2);
651 trans_elastic.commonData.spatialPositions = "DISPLACEMENT";
652 trans_elastic.commonData.meshPositions = "MESH_NODE_POSITIONS";
653 std::map<int,boost::shared_ptr<SmallStrainTranverslyIsotropicADouble> > tranversly_isotropic_adouble_ptr_map;
654 std::map<int,boost::shared_ptr<SmallStrainTranverslyIsotropicDouble> > tranversly_isotropic_double_ptr_map;
655 bool trans_iso_blocks = false;
656 for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
657 //Get block name
658 string name = it->getName();
659 if (name.compare(0,20,"MAT_ELASTIC_TRANSISO") == 0) {
660 cout<<"================================ it is MAT_ELASTIC_TRANSISO "<<endl;
661 trans_iso_blocks = true;
662 int id = it->getMeshsetId();
663 Mat_Elastic_TransIso mydata;
664 ierr = it->getAttributeDataStructure(mydata); CHKERRQ(ierr);
665 tranversly_isotropic_adouble_ptr_map[id] = boost::make_shared<SmallStrainTranverslyIsotropicADouble>();
666 tranversly_isotropic_double_ptr_map[id] = boost::make_shared<SmallStrainTranverslyIsotropicDouble>();
667 //nu_p, nu_pz, E_p, E_z, G_zp
668 tranversly_isotropic_adouble_ptr_map.at(id)->E_p = mydata.data.Youngp;
669 tranversly_isotropic_double_ptr_map.at(id)->E_p = mydata.data.Youngp;
670 tranversly_isotropic_adouble_ptr_map.at(id)->E_z = mydata.data.Youngz;
671 tranversly_isotropic_double_ptr_map.at(id)->E_z = mydata.data.Youngz;
672 tranversly_isotropic_adouble_ptr_map.at(id)->nu_p = mydata.data.Poissonp;
673 tranversly_isotropic_double_ptr_map.at(id)->nu_p = mydata.data.Poissonp;
674 tranversly_isotropic_adouble_ptr_map.at(id)->nu_pz = mydata.data.Poissonpz;
675 tranversly_isotropic_double_ptr_map.at(id)->nu_pz = mydata.data.Poissonpz;
676 double shear_zp;
677 if(mydata.data.Shearzp!=0) {
678 shear_zp = mydata.data.Shearzp;
679 } else {
680 shear_zp = mydata.data.Youngz/(2*(1+mydata.data.Poissonpz));
681 }
682 tranversly_isotropic_adouble_ptr_map.at(it->getMeshsetId())->G_zp = shear_zp;
683 tranversly_isotropic_double_ptr_map.at(it->getMeshsetId())->G_zp = shear_zp;
684 //get tets from block where material is defined
685 EntityHandle meshset = it->getMeshset();
686 rval = m_field.get_moab().get_entities_by_type(
687 meshset,MBTET,trans_elastic.setOfBlocks[id].tEts,true
688 ); CHKERRQ_MOAB(rval);
689
690// cout<<"================================ trans_elastic.setOfBlocks[id].tEts "<<trans_elastic.setOfBlocks[id].tEts.size()<<endl;
691 //intersection of new and old tets for tran-iso part
692 trans_elastic.setOfBlocks[id].tEts = intersect(trans_elastic.setOfBlocks[id].tEts,LatestRefinedTets);
693// cout<<"================================ trans_elastic.setOfBlocks[id].tEts "<<trans_elastic.setOfBlocks[id].tEts.size()<<endl;
694
695 //adding material to nonlinear class
696 trans_elastic.setOfBlocks[id].iD = id;
697 //note that material parameters are defined internally in material model
698 trans_elastic.setOfBlocks[id].E = 0; // this is not working for this material
699 trans_elastic.setOfBlocks[id].PoissonRatio = 0; // this is not working for this material
700 trans_elastic.setOfBlocks[id].materialDoublePtr = tranversly_isotropic_double_ptr_map.at(id);
701 trans_elastic.setOfBlocks[id].materialAdoublePtr = tranversly_isotropic_adouble_ptr_map.at(id);
702 ierr = m_field.seed_finite_elements(trans_elastic.setOfBlocks[id].tEts); CHKERRQ(ierr);
703 }
704 }
705
706
707 ierr = m_field.add_finite_element("TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
708 ierr = m_field.add_finite_element("TRAN_ISOTROPIC_ELASTIC",MF_ZERO); CHKERRQ(ierr);
709 ierr = m_field.modify_finite_element_add_field_row("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
710 ierr = m_field.modify_finite_element_add_field_col("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
711 ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
712 ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","POTENTIAL_FIELD"); CHKERRQ(ierr);
713
714 if(m_field.check_field("MESH_NODE_POSITIONS")) {
715 ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
716 }
717 for(
718 map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
719 sit!=trans_elastic.setOfBlocks.end();sit++
720 ) {
721
722
723 cout<<" sit->second.tEts "<<sit->second.tEts.size()<<endl;
724 ierr = m_field.add_ents_to_finite_element_by_type(sit->second.tEts,MBTET,"TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
725 }
726
727
728 //Rhs
729 trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("DISPLACEMENT",trans_elastic.commonData));
730 trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("POTENTIAL_FIELD",trans_elastic.commonData));
731 if(m_field.check_field("MESH_NODE_POSITIONS")) {
732 trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("MESH_NODE_POSITIONS",trans_elastic.commonData));
733 }
734 map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
735 for(;sit!=trans_elastic.setOfBlocks.end();sit++) {
736 trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress("DISPLACEMENT",sit->second,trans_elastic.commonData,2,false,false,true));
737 trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpRhsPiolaKirchhoff("DISPLACEMENT",sit->second,trans_elastic.commonData));
738 }
739
740 //Lhs
741 trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("DISPLACEMENT",trans_elastic.commonData));
742 trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("POTENTIAL_FIELD",trans_elastic.commonData));
743 if(m_field.check_field("MESH_NODE_POSITIONS")) {
744 trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("MESH_NODE_POSITIONS",trans_elastic.commonData));
745 }
746 sit = trans_elastic.setOfBlocks.begin();
747 for(;sit!=trans_elastic.setOfBlocks.end();sit++) {
748 trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress("DISPLACEMENT",sit->second,trans_elastic.commonData,2,true,false,true));
749 trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpLhsPiolaKirchhoff_dx("DISPLACEMENT","DISPLACEMENT",sit->second,trans_elastic.commonData));
750 }
751
752
753
754
755 //================================================================================================================================
756 // Interface Cohesive elemetns
757 //================================================================================================================================
758 //FE Interface
759
760 ierr = m_field.add_finite_element("INTERFACE"); CHKERRQ(ierr);
761 ierr = m_field.modify_finite_element_add_field_row("INTERFACE","DISPLACEMENT"); CHKERRQ(ierr);
762 ierr = m_field.modify_finite_element_add_field_col("INTERFACE","DISPLACEMENT"); CHKERRQ(ierr);
763 ierr = m_field.modify_finite_element_add_field_data("INTERFACE","DISPLACEMENT"); CHKERRQ(ierr);
764 ierr = m_field.modify_finite_element_add_field_data("INTERFACE","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
765
766 boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation> interface_materials;
767
768 //FIXME this in fact allow only for one type of interface,
769 //problem is Young Moduls in interface mayoung_modulusterial
770
771
772 double interface_beta, interface_ft, interface_Gf, interface_h;
773 ierr = PetscOptionsGetReal(0,"-interface_beta",&interface_beta,0); CHKERRQ(ierr);
774 ierr = PetscOptionsGetReal(0,"-interface_ft",&interface_ft,0); CHKERRQ(ierr);
775 ierr = PetscOptionsGetReal(0,"-interface_Gf",&interface_Gf,0); CHKERRQ(ierr);
776 ierr = PetscOptionsGetReal(0,"-interface_h",&interface_h,0); CHKERRQ(ierr);
777
778
779 for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
780 cout << endl << *it << endl;
781 //Get block name
782 string name = it->getName();
783 if (name.compare(0,10,"MAT_INTERF") == 0) {
784 Mat_Interf mydata;
785 ierr = it->getAttributeDataStructure(mydata); CHKERRQ(ierr);
786 cout << mydata;
787 interface_materials.push_back(new CohesiveInterfaceElement::PhysicalEquation(m_field));
788// interface_materials.back().h = 0.134; // we used this because E0=35Gpa and Em=4.7Gpa
789 interface_materials.back().youngModulus = mydata.data.alpha;
790// interface_materials.back().beta = mydata.data.beta;
791// interface_materials.back().ft = mydata.data.ft;
792// interface_materials.back().Gf = mydata.data.Gf;
793
794 interface_materials.back().h = interface_h;
795 interface_materials.back().beta = interface_beta;
796 interface_materials.back().ft = interface_ft;
797 interface_materials.back().Gf = interface_Gf;
798
799 EntityHandle meshset = it->getMeshset();
800 Range tris; //All the triangles after splitting
801 rval = moab.get_entities_by_type(meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
802// cout<<"tris.size() "<< tris.size() <<endl;
803 Range ents3d;
804 rval = moab.get_adjacencies(tris,3,false,ents3d,moab::Interface::UNION); CHKERRQ_MOAB(rval);
805 interface_materials.back().pRisms = ents3d.subset_by_type(MBPRISM);
806// cout<<"unite( interface_prisms, ents3d.subset_by_type(MBPRISM) ) "<<unite( interface_prisms, ents3d.subset_by_type(MBPRISM) )<<endl;
807 }
808 }
809
810 ierr = m_field.seed_finite_elements(interface_materials.back().pRisms); CHKERRQ(ierr);
811// cout<<"interface_materials.back().pRisms.size() "<< interface_materials.back().pRisms <<endl;
812 ierr = m_field.add_ents_to_finite_element_by_type(interface_materials.back().pRisms,MBPRISM,"INTERFACE"); CHKERRQ(ierr);
813
814
815 // cout<<"young_modulus "<<young_modulus<<endl;
816 { //FIXME
817 boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>::iterator pit = interface_materials.begin();
818 for(; pit != interface_materials.end();pit++) {
819 pit->youngModulus = young_modulus; //Young's modulus of bulk elastic material (then for interface E0=E/h)
820 }
821 }
822
823 CohesiveInterfaceElement cohesive_elements(m_field); //make it golbal, so that its other uses can see it
824 ierr = cohesive_elements.addOps("DISPLACEMENT",interface_materials); CHKERRQ(ierr);
825 //================================================================================================================================
826
827
828 BCs_RVELagrange_Periodic lagrangian_element_periodic(m_field);
829// BCs_RVELagrange_Trac lagrangian_element_trac(m_field);
830
831 lagrangian_element_periodic.addLagrangiangElement("LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",preriodic_prisms);
832
833
834 //Adding element LAGRANGE_ELEM_TRANS
835// lagrangian_element_trac.addLagrangiangElement("LAGRANGE_ELEM_TRANS","DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS","MESH_NODE_POSITIONS");
836 ierr = m_field.add_finite_element("LAGRANGE_ELEM_TRANS",MF_ZERO); CHKERRQ(ierr);
837 //============================================================================================================
838 //C row as Lagrange_mul_disp and col as DISPLACEMENT
839 ierr = m_field.modify_finite_element_add_field_row("LAGRANGE_ELEM_TRANS","LAGRANGE_MUL_RIGID_TRANS"); CHKERRQ(ierr);
840 ierr = m_field.modify_finite_element_add_field_col("LAGRANGE_ELEM_TRANS","DISPLACEMENT"); CHKERRQ(ierr);
841 //CT col as Lagrange_mul_disp and row as DISPLACEMENT
842 ierr = m_field.modify_finite_element_add_field_col("LAGRANGE_ELEM_TRANS","LAGRANGE_MUL_RIGID_TRANS"); CHKERRQ(ierr);
843 ierr = m_field.modify_finite_element_add_field_row("LAGRANGE_ELEM_TRANS","DISPLACEMENT"); CHKERRQ(ierr);
844 //data
845 ierr = m_field.modify_finite_element_add_field_data("LAGRANGE_ELEM_TRANS","LAGRANGE_MUL_RIGID_TRANS"); CHKERRQ(ierr);
846 ierr = m_field.modify_finite_element_add_field_data("LAGRANGE_ELEM_TRANS","DISPLACEMENT"); CHKERRQ(ierr);
847 //============================================================================================================
848
849
850 Range Tri103;
851 for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
852 if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
853
854 //total tri including old + new
855 rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,Tri103,true); CHKERRQ_MOAB(rval);
856 ierr = PetscPrintf(PETSC_COMM_WORLD,"Old triangles on SideSet 103 = %d\n",Tri103.size()); CHKERRQ(ierr);
857 //Only triangles in the new bit-level
858 Tri103 = intersect(AllTris, Tri103);
859 ierr = PetscPrintf(PETSC_COMM_WORLD,"New triangles on SideSet 103 = %d\n",Tri103.size()); CHKERRQ(ierr);
860 ierr = m_field.add_ents_to_finite_element_by_type(Tri103,MBTRI,"LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
861 }
862 }
863
864// build finite elements
865 ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
866 //build adjacencies
867 ierr = m_field.build_adjacencies(problem_bit_level); CHKERRQ(ierr);
868
869
870 DMType dm_name = "PLASTIC_PROB";
871 ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
872
873 DM dm;
874 ierr = DMCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(ierr);
875 ierr = DMSetType(dm,dm_name);CHKERRQ(ierr);
876
877 //set dm datastruture which created mofem datastructures
878 ierr = DMMoFEMCreateMoFEM(dm,&m_field,dm_name,problem_bit_level); CHKERRQ(ierr);
879 ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
880 ierr = DMMoFEMSetIsPartitioned(dm,is_partitioned); CHKERRQ(ierr);
881
882 //add elements to dm
883 ierr = DMMoFEMAddElement(dm,"PLASTIC"); CHKERRQ(ierr);
884 ierr = DMMoFEMAddElement(dm,"TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
885 ierr = DMMoFEMAddElement(dm,"INTERFACE"); CHKERRQ(ierr);
886 ierr = DMMoFEMAddElement(dm,"LAGRANGE_ELEM"); CHKERRQ(ierr);
887 ierr = DMMoFEMAddElement(dm,"LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
888 ierr = DMSetUp(dm); CHKERRQ(ierr);
889
890 //create matrices
891 Vec F,D;
892 Mat Aij;
893 ierr = DMCreateGlobalVector_MoFEM(dm,&D); CHKERRQ(ierr);
894 ierr = VecDuplicate(D,&F); CHKERRQ(ierr);
895 ierr = DMCreateMatrix_MoFEM(dm,&Aij); CHKERRQ(ierr);
896 ierr = VecZeroEntries(D); CHKERRQ(ierr);
897 ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
898 ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
899 ierr = DMoFEMMeshToLocalVector(dm,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
900 ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
901
902 vector<Vec> Fvec(6); //jthis will be used to caluclate the homogenised stiffness matix
903 for(int ii = 0;ii<6;ii++) {
904 ierr = VecDuplicate(D,&Fvec[ii]); CHKERRQ(ierr);
905 ierr = VecZeroEntries(Fvec[ii]); CHKERRQ(ierr);
906 ierr = VecGhostUpdateBegin(Fvec[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
907 ierr = VecGhostUpdateEnd(Fvec[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
908 }
909
910 //assemble Aij and F
911 SmallStrainPlasticity small_strain_plasticity(m_field);
912 {
913 PetscBool bbar = PETSC_TRUE;
914 ierr = PetscOptionsGetBool(0,"-my_bbar",&bbar,0); CHKERRQ(ierr);
915 small_strain_plasticity.commonData.bBar = bbar;
916 }
917 {
918 small_strain_plasticity.feRhs.getOpPtrVector().push_back(
919 new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
920 "DISPLACEMENT",small_strain_plasticity.commonData
921 )
922 );
923 small_strain_plasticity.feRhs.getOpPtrVector().push_back(
924 new SmallStrainPlasticity::OpCalculateStress(
925 m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp,false
926 )
927 );
928 small_strain_plasticity.feRhs.getOpPtrVector().push_back(
929 new SmallStrainPlasticity::OpAssembleRhs(
930 "DISPLACEMENT",small_strain_plasticity.commonData
931 )
932 );
933 small_strain_plasticity.feLhs.getOpPtrVector().push_back(
934 new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
935 "DISPLACEMENT",small_strain_plasticity.commonData
936 )
937 );
938 small_strain_plasticity.feLhs.getOpPtrVector().push_back(
939 new SmallStrainPlasticity::OpCalculateStress(
940 m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp,false
941 )
942 );
943 small_strain_plasticity.feLhs.getOpPtrVector().push_back(
944 new SmallStrainPlasticity::OpAssembleLhs(
945 "DISPLACEMENT",small_strain_plasticity.commonData
946 )
947 );
948 small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
949 new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
950 "DISPLACEMENT",small_strain_plasticity.commonData
951 )
952 );
953 small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
954 new SmallStrainPlasticity::OpCalculateStress(
955 m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp,false
956 )
957 );
958 small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
959 new SmallStrainPlasticity::OpUpdate(
960 "DISPLACEMENT",small_strain_plasticity.commonData
961 )
962 );
963 ierr = small_strain_plasticity.createInternalVariablesTag(); CHKERRQ(ierr);
964 }
965
966
967
968 lagrangian_element_periodic.setRVEBCsOperatorsNonlinear("DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",Aij,Fvec,F,given_strain);
969 BCs_RVELagrange_Trac_Rigid_Trans lagrangian_element_rigid_body_trans(m_field);
970 lagrangian_element_rigid_body_trans.setRVEBCsRigidBodyTranOperators("DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS",Aij, lagrangian_element_periodic.setOfRVEBC);
971
972
973 TimeForceScale time_force_scale("-my_load_history",false);
974 lagrangian_element_periodic.methodsOp.push_back(new TimeForceScale("-my_macro_strian_history",false));
975
976 //Adding elements to DMSnes
977 //Rhs
978 ierr = DMMoFEMSNESSetFunction(dm,"PLASTIC",&small_strain_plasticity.feRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
979 ierr = DMMoFEMSNESSetFunction(dm,"TRAN_ISOTROPIC_ELASTIC",&trans_elastic.feRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
980 ierr = DMMoFEMSNESSetFunction(dm,"INTERFACE",&cohesive_elements.feRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
981 ierr = DMMoFEMSNESSetFunction(dm,"LAGRANGE_ELEM",&lagrangian_element_periodic.feRVEBCRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
982 ierr = DMMoFEMSNESSetFunction(dm,"LAGRANGE_ELEM",&lagrangian_element_periodic.feRVEBCRhsResidual,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
983
984 //Lhs
985 ierr = DMMoFEMSNESSetJacobian(dm,"PLASTIC",&small_strain_plasticity.feLhs,NULL,NULL); CHKERRQ(ierr);
986 ierr = DMMoFEMSNESSetJacobian(dm,"TRAN_ISOTROPIC_ELASTIC",&trans_elastic.feLhs,NULL,NULL); CHKERRQ(ierr);
987 ierr = DMMoFEMSNESSetJacobian(dm,"INTERFACE",&cohesive_elements.feLhs,NULL,NULL); CHKERRQ(ierr);
988 ierr = DMMoFEMSNESSetJacobian(dm,"LAGRANGE_ELEM",&lagrangian_element_periodic.feRVEBCLhs,NULL,NULL); CHKERRQ(ierr);
989 ierr = DMMoFEMSNESSetJacobian(dm,"LAGRANGE_ELEM_TRANS",&lagrangian_element_rigid_body_trans.feRVEBCLhs,NULL,NULL); CHKERRQ(ierr);
990
991 // Create Time Stepping solver
992 SNES snes;
993 SnesCtx *snes_ctx;
994 SNES snes_one_step_only;
995
996
997 ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
998 //ierr = SNESSetDM(snes,dm); CHKERRQ(ierr);
999 ierr = DMMoFEMGetSnesCtx(dm,&snes_ctx); CHKERRQ(ierr);
1000 ierr = SNESSetFunction(snes,F,SnesRhs,snes_ctx); CHKERRQ(ierr);
1001 ierr = SNESSetJacobian(snes,Aij,Aij,SnesMat,snes_ctx); CHKERRQ(ierr);
1002 ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
1003
1004 ierr = SNESCreate(PETSC_COMM_WORLD,&snes_one_step_only); CHKERRQ(ierr);
1005 ierr = SNESSetFunction(snes_one_step_only,F,SnesRhs,snes_ctx); CHKERRQ(ierr);
1006 ierr = SNESSetJacobian(snes_one_step_only,Aij,Aij,SnesMat,snes_ctx); CHKERRQ(ierr);
1007 ierr = SNESSetFromOptions(snes_one_step_only); CHKERRQ(ierr);
1008 double atol,rtol,stol;
1009 int maxit,maxf;
1010 SNESGetTolerances(snes_one_step_only,&atol,&rtol,&stol,&maxit,&maxf);
1011 maxit = 1;
1012 atol *= 0;
1013 SNESSetTolerances(snes_one_step_only,atol,rtol,stol,maxit,maxf);
1014 SNESLineSearch linesearch;
1015 SNESGetLineSearch(snes_one_step_only,&linesearch);
1016 SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC);
1017
1018
1019
1020 PostProcVolumeOnRefinedMesh post_proc(m_field);
1021 ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
1022 ierr = post_proc.addFieldValuesPostProc("DISPLACEMENT"); CHKERRQ(ierr);
1023 ierr = post_proc.addFieldValuesGradientPostProc("DISPLACEMENT"); CHKERRQ(ierr);
1024 ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
1025
1026 // Volume calculation
1027 //==============================================================================================================================
1028 VolumeElementForcesAndSourcesCore vol_elem(m_field);
1029 Vec volume_vec;
1030 int volume_vec_ghost[] = { 0 };
1031 ierr = VecCreateGhost(
1032 PETSC_COMM_WORLD,(!m_field.get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
1033 ); CHKERRQ(ierr);
1034 ierr = VecZeroEntries(volume_vec); CHKERRQ(ierr);
1035 vol_elem.getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
1036
1037 ierr = m_field.loop_finite_elements("PLASTIC_PROB","PLASTIC", vol_elem); CHKERRQ(ierr);
1038 ierr = m_field.loop_finite_elements("PLASTIC_PROB","TRAN_ISOTROPIC_ELASTIC", vol_elem); CHKERRQ(ierr);
1039
1040 ierr = VecAssemblyBegin(volume_vec); CHKERRQ(ierr);
1041 ierr = VecAssemblyEnd(volume_vec); CHKERRQ(ierr);
1042 double rve_volume;
1043 ierr = VecSum(volume_vec,&rve_volume); CHKERRQ(ierr);
1044
1045 // ierr = VecView(volume_vec,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
1046 ierr = PetscPrintf(PETSC_COMM_WORLD,"RVE Volume %3.2g\n",rve_volume); CHKERRQ(ierr);
1047 ierr = VecDestroy(&volume_vec);
1048 //=============================================================================================================================
1049
1050 double final_time = 1,delta_time = 0.1;
1051 ierr = PetscOptionsGetReal(0,"-my_final_time",&final_time,0); CHKERRQ(ierr);
1052 ierr = PetscOptionsGetReal(0,"-my_delta_time",&delta_time,0); CHKERRQ(ierr);
1053 double delta_time0 = delta_time;
1054
1055// ierr = MatView(Aij,PETSC_VIEWER_DRAW_SELF); CHKERRQ(ierr);
1056// std::string wait;
1057// std::cin >> wait;
1058
1059 Vec D0;
1060 ierr = VecDuplicate(D,&D0); CHKERRQ(ierr);
1061
1062 int step = 0;
1063 double t = 0;
1064 SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
1065 for(;t<final_time;step++) {
1066 t += delta_time;
1067
1068// if(t>final_time){
1069// break;
1070// }
1071 PetscPrintf(
1072 PETSC_COMM_WORLD,"Step %d Time %6.4g final time %3.2g\n",step,t,final_time
1073 );
1074 //set time
1075 lagrangian_element_periodic.getLoopFeRVEBCRhs().ts_t = t;
1076 //solve problem at step
1077 ierr = VecAssemblyBegin(D);
1078 ierr = VecAssemblyEnd(D);
1079 ierr = VecCopy(D,D0); CHKERRQ(ierr);
1080 if(step == 0 || reason < 0) {
1081 ierr = SNESSetLagJacobian(snes,-2); CHKERRQ(ierr);
1082 } else {
1083 ierr = SNESSetLagJacobian(snes,-1); CHKERRQ(ierr);
1084 }
1085 ierr = SNESSetLagJacobian(snes_one_step_only,-2);
1086 ierr = SNESSolve(snes_one_step_only,PETSC_NULL,D); CHKERRQ(ierr);
1087 ierr = SNESGetConvergedReason(snes_one_step_only,&reason); CHKERRQ(ierr);
1088 if(reason<0) {
1089 ierr = SNESSolve(snes,PETSC_NULL,D); CHKERRQ(ierr);
1090 }
1091 int its;
1092 ierr = SNESGetIterationNumber(snes,&its); CHKERRQ(ierr);
1093
1094 ierr = PetscPrintf(
1095 PETSC_COMM_WORLD,"number of Newton iterations = %D\n",its
1096 ); CHKERRQ(ierr);
1097
1098 ierr = SNESGetConvergedReason(snes,&reason); CHKERRQ(ierr);
1099
1100 if(reason<0) {
1101 t -= delta_time;
1102 delta_time *= 0.1;
1103 ierr = VecCopy(D0,D); CHKERRQ(ierr);
1104 } else {const int its_d = 60;
1105 const double gamma = 0.5;
1106 const double max_reudction = 1;
1107 const double min_reduction = 1e-1;
1108 double reduction;
1109 reduction = pow((double)its_d/(double)(its+1),gamma);
1110 if(delta_time >= max_reudction*delta_time0 && reduction > 1) {
1111 reduction = 1;
1112 } else if(delta_time <= min_reduction*delta_time0 && reduction < 1) {
1113 reduction = 1;
1114 }
1115
1116 ierr = PetscPrintf(
1117 PETSC_COMM_WORLD,
1118 "reduction delta_time = %6.4e delta_time = %6.4e\n",
1119 reduction,delta_time
1120 ); CHKERRQ(ierr);
1121 delta_time *= reduction;
1122 if(reduction>1 && delta_time < min_reduction*delta_time0) {
1123 delta_time = min_reduction*delta_time0;
1124 }
1125
1126 ierr = DMoFEMMeshToGlobalVector(
1127 dm,D,INSERT_VALUES,SCATTER_REVERSE
1128 ); CHKERRQ(ierr);
1129
1130 ierr = DMoFEMLoopFiniteElements(
1131 dm,"PLASTIC",&small_strain_plasticity.feUpdate
1132 ); CHKERRQ(ierr);
1133
1134
1135 ierr = DMoFEMLoopFiniteElements(
1136 dm,"INTERFACE",&cohesive_elements.feHistory
1137 ); CHKERRQ(ierr);
1138
1139 //Save data on mesh
1140 {
1141 ierr = DMoFEMLoopFiniteElements(
1142 dm,"PLASTIC",&post_proc
1143 ); CHKERRQ(ierr);
1144
1145 string out_file_name;
1146 {
1147 std::ostringstream stm;
1148 stm << "out_matrix_" << step << ".h5m";
1149 out_file_name = stm.str();
1150 }
1151
1152 ierr = PetscPrintf(
1153 PETSC_COMM_WORLD,"out file %s\n",out_file_name.c_str()
1154 ); CHKERRQ(ierr);
1155
1156 rval = post_proc.postProcMesh.write_file(
1157 out_file_name.c_str(),"MOAB","PARALLEL=WRITE_PART"
1158 ); CHKERRQ_MOAB(rval);
1159
1160
1161
1162 ierr = DMoFEMLoopFiniteElements(
1163 dm,"TRAN_ISOTROPIC_ELASTIC",&post_proc
1164 ); CHKERRQ(ierr);
1165 {
1166 std::ostringstream stm;
1167 stm << "out_fibres_" << step << ".h5m";
1168 out_file_name = stm.str();
1169 }
1170 ierr = PetscPrintf(
1171 PETSC_COMM_WORLD,"out file %s\n",out_file_name.c_str()
1172 ); CHKERRQ(ierr);
1173
1174 rval = post_proc.postProcMesh.write_file(
1175 out_file_name.c_str(),"MOAB","PARALLEL=WRITE_PART"
1176 ); CHKERRQ_MOAB(rval);
1177
1178 }
1179
1180 //===================================== Homgenised stress (for given strain) ====================================================
1181 VectorDouble StressHomo;
1182 StressHomo.resize(6);
1183 StressHomo.clear();
1184
1185 //calculate honmogenised stress
1186 //create a vector for 6 components of homogenized stress
1187 Vec stress_homo;
1188 int stress_homo_ghost[] = { 0,1,2,3,4,5,6 };
1189 ierr = VecCreateGhost(
1190 PETSC_COMM_WORLD,(!m_field.get_comm_rank())?6:0,6,6,stress_homo_ghost,&stress_homo
1191 ); CHKERRQ(ierr);
1192
1193 lagrangian_element_periodic.setRVEBCsHomoStressOperators("DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",stress_homo);
1194
1195
1196 PetscScalar *avec;
1197 ierr = VecGetArray(stress_homo,&avec); CHKERRQ(ierr);
1198 ierr = VecZeroEntries(stress_homo); CHKERRQ(ierr);
1199 ierr = m_field.loop_finite_elements(
1200 "PLASTIC_PROB","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCStress()
1201 ); CHKERRQ(ierr);
1202 ierr = PetscOptionsGetReal(
1203 PETSC_NULL,"-my_rve_volume",&rve_volume,PETSC_NULL
1204 ); CHKERRQ(ierr);
1205 ierr = VecAssemblyBegin(stress_homo); CHKERRQ(ierr);
1206 ierr = VecAssemblyEnd(stress_homo); CHKERRQ(ierr);
1207 ierr = VecScale(stress_homo,1.0/rve_volume); CHKERRQ(ierr);
1208 ierr = VecGhostUpdateBegin(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1209 ierr = VecGhostUpdateEnd(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1210
1211 for(int jj=0; jj<6; jj++) {
1212 StressHomo(jj) = avec[jj];
1213 }
1214
1215 double scale;
1216 ierr = time_force_scale.getForceScale(t,scale); CHKERRQ(ierr);
1217
1218 PetscPrintf(PETSC_COMM_WORLD,
1219 "Macro_Strain Homo_Stress Path "
1220 );
1221
1222 //cout<<"Macro-strain ";
1223 for(int ii=0; ii<6; ii++) {
1224 PetscPrintf(
1225 PETSC_COMM_WORLD,
1226 "%8.5e\t",
1227 t*given_strain(ii)
1228 );
1229 }
1230
1231 //cout<<"Homo stress ";
1232 for(int ii=0; ii<6; ii++) {
1233 PetscPrintf(
1234 PETSC_COMM_WORLD,
1235 "%8.5e\t",
1236 StressHomo(ii)
1237 );
1238 }
1239
1240 PetscPrintf(PETSC_COMM_WORLD,
1241 "\n"
1242 );
1243// //==============================================================================================================================
1244 }
1245 }
1246
1247 ierr = VecDestroy(&D0); CHKERRQ(ierr);
1248 ierr = MatDestroy(&Aij); CHKERRQ(ierr);
1249 ierr = VecDestroy(&F); CHKERRQ(ierr);
1250 ierr = VecDestroy(&D); CHKERRQ(ierr);
1251
1252 }
1253 CATCH_ERRORS;
1254
1255
1256 MoFEM::Core::Finalize();
1257}
BCs_RVELagrange_Disp.hpp
BCs_RVELagrange_Periodic.hpp
BCs_RVELagrange_Trac.hpp
BCs_RVELagrange_Trac_Rigid_Trans.hpp
CohesiveInterfaceElement.hpp
Implementation of linear interface element.
NonLinearElasticElement.hpp
Operators and data structures for non-linear elastic analysis.
PostProcOnRefMesh.hpp
Post-process fields on refined mesh.
SmallStrainPlasticity.hpp
Operators and data structures for small strain plasticity.
SmallStrainPlasticityMaterialModels.hpp
Small Strain plasticity material models.
VolumeCalculation.hpp
Operator can be used with any volume element to calculate sum of volumes of all volumes in the set.
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
NOBASE
@ NOBASE
Definition: definitions.h:72
MOAB_THROW
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:554
NOFIELD
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:97
H1
@ H1
continuous field
Definition: definitions.h:98
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
SIDESET
@ SIDESET
Definition: definitions.h:160
UNKNOWNNAME
@ UNKNOWNNAME
Definition: definitions.h:171
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:467
MU
#define MU(E, NU)
Definition: fem_tools.h:33
LAMBDA
#define LAMBDA(E, NU)
Definition: fem_tools.h:32
poisson_ratio
auto poisson_ratio
Definition: fluid_structure_eigenproblem.cpp:71
order
int order
Definition: fluid_structure_eigenproblem.cpp:68
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1082
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: DMMMoFEM.cpp:677
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: DMMMoFEM.cpp:130
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:462
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: DMMMoFEM.cpp:482
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: DMMMoFEM.cpp:718
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMMoFEM.cpp:1156
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
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: DMMMoFEM.cpp:545
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMMoFEM.cpp:1126
MoFEM::DMMoFEMGetSnesCtx
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMMoFEM.cpp:1053
MoFEM::DMoFEMMeshToGlobalVector
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:494
MoFEM::CoreInterface::modify_finite_element_add_field_row
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
MoFEM::CoreInterface::modify_finite_element_add_field_col
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::set_field_order
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
MoFEM::CoreInterface::add_ents_to_field_by_type
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
PostProcTemplateOnRefineMesh::addFieldValuesGradientPostProc
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
Definition: PostProcOnRefMesh.hpp:207
PostProcTemplateOnRefineMesh::addFieldValuesPostProc
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
Definition: PostProcOnRefMesh.hpp:165
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.
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:59
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:81
Face_CenPos_Handle_multiIndex
multi_index_container< Face_CenPos_Handle, indexed_by< ordered_non_unique< tag< xcoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::xcoord > >, ordered_non_unique< tag< ycoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::ycoord > >, ordered_non_unique< tag< zcoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::zcoord > >, ordered_unique< tag< Tri_Hand_tag >, member< Face_CenPos_Handle, const EntityHandle,&Face_CenPos_Handle::Tri_Hand > >, ordered_unique< tag< Composite_xyzcoord >, composite_key< Face_CenPos_Handle, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::xcoord >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::ycoord >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::zcoord > > > > > Face_CenPos_Handle_multiIndex
Definition: homogenisation_periodic_atom_test.cpp:96
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:35
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:78
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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:148
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:153
MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:163
MoFEM::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:39
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:193
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1096
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:183
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:32
std
STL namespace.
young_modulus
double young_modulus
Definition: plastic.cpp:105
scale
double scale
Definition: plastic.cpp:103
t
constexpr double t
plate stiffness
Definition: plate.cpp:76
Face_CenPos_Handle_multiIndex
multi_index_container< Face_CenPos_Handle, indexed_by< ordered_non_unique< tag< xcoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::xcoord > >, ordered_non_unique< tag< ycoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::ycoord > >, ordered_non_unique< tag< zcoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::zcoord > >, ordered_unique< tag< Tri_Hand_tag >, member< Face_CenPos_Handle, const EntityHandle,&Face_CenPos_Handle::Tri_Hand > >, ordered_unique< tag< Composite_xyzcoord >, composite_key< Face_CenPos_Handle, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::xcoord >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::ycoord >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::zcoord > > > > > Face_CenPos_Handle_multiIndex
Definition: rve_mech_plas_interface_periodic.cpp:115
main
int main(int argc, char *argv[])
Definition: rve_mech_plas_interface_periodic.cpp:122
help
static char help[]
Definition: rve_mech_plas_interface_periodic.cpp:71
BCs_RVELagrange_Disp::feRVEBCLhs
MyTriFE feRVEBCLhs
Definition: BCs_RVELagrange_Disp.hpp:30
BCs_RVELagrange_Disp::setOfRVEBC
map< int, RVEBC_Data > setOfRVEBC
maps side set id with appropriate FluxData
Definition: BCs_RVELagrange_Disp.hpp:56
BCs_RVELagrange_Periodic
Definition: BCs_RVELagrange_Periodic.hpp:19
BCs_RVELagrange_Periodic::getLoopFeRVEBCRhs
MyPrismFE & getLoopFeRVEBCRhs()
Definition: BCs_RVELagrange_Periodic.hpp:39
BCs_RVELagrange_Periodic::feRVEBCRhs
MyPrismFE feRVEBCRhs
Definition: BCs_RVELagrange_Periodic.hpp:33
BCs_RVELagrange_Periodic::methodsOp
boost::ptr_vector< MethodForForceScaling > methodsOp
Definition: BCs_RVELagrange_Periodic.hpp:21
BCs_RVELagrange_Periodic::setRVEBCsHomoStressOperators
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec stress_homo)
Definition: BCs_RVELagrange_Periodic.hpp:1312
BCs_RVELagrange_Periodic::getLoopFeRVEBCStress
MyPrismFE & getLoopFeRVEBCStress()
Definition: BCs_RVELagrange_Periodic.hpp:42
BCs_RVELagrange_Periodic::feRVEBCLhs
MyPrismFE feRVEBCLhs
Definition: BCs_RVELagrange_Periodic.hpp:34
BCs_RVELagrange_Periodic::addLagrangiangElement
PetscErrorCode addLagrangiangElement(const string element_name, const string field_name, const string lagrang_field_name, const string mesh_nodals_positions, Range &periodic_prisms)
Definition: BCs_RVELagrange_Periodic.hpp:52
BCs_RVELagrange_Periodic::feRVEBCRhsResidual
MyPrismFE feRVEBCRhsResidual
Definition: BCs_RVELagrange_Periodic.hpp:36
BCs_RVELagrange_Periodic::setRVEBCsOperatorsNonlinear
PetscErrorCode setRVEBCsOperatorsNonlinear(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &fvec, Vec f, VectorDouble given_strain)
Definition: BCs_RVELagrange_Periodic.hpp:1091
BCs_RVELagrange_Trac_Rigid_Trans
Definition: BCs_RVELagrange_Trac_Rigid_Trans.hpp:18
BCs_RVELagrange_Trac_Rigid_Trans::setRVEBCsRigidBodyTranOperators
PetscErrorCode setRVEBCsRigidBodyTranOperators(string field_name, string lagrang_field_name, Mat aij, map< int, RVEBC_Data > setOfRVEBC)
Definition: BCs_RVELagrange_Trac_Rigid_Trans.hpp:109
CohesiveElement::CohesiveInterfaceElement::addOps
MoFEMErrorCode addOps(const std::string field_name, boost::ptr_vector< CohesiveInterfaceElement::PhysicalEquation > &interfaces)
Driver function settting all operators needed for interface element.
Definition: CohesiveInterfaceElement.hpp:580
Composite_xyzcoord
Definition: homogenisation_periodic_atom_test.cpp:72
Face_CenPos_Handle
Definition: homogenisation_periodic_atom_test.cpp:61
Face_CenPos_Handle::zcoord
double zcoord
Definition: homogenisation_periodic_atom_test.cpp:62
Face_CenPos_Handle::Tri_Hand
const EntityHandle Tri_Hand
Definition: homogenisation_periodic_atom_test.cpp:63
Face_CenPos_Handle::Face_CenPos_Handle
Face_CenPos_Handle(double _xcoord, double _ycoord, double _zcoord, const EntityHandle _Tri_Hand)
Definition: rve_mech_plas_interface_periodic.cpp:83
Face_CenPos_Handle::xcoord
double xcoord
Definition: homogenisation_periodic_atom_test.cpp:62
Face_CenPos_Handle::ycoord
double ycoord
Definition: homogenisation_periodic_atom_test.cpp:62
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
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:92
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:26
MoFEM::DeprecatedCoreInterface::seed_ref_level
DEPRECATED MoFEMErrorCode seed_ref_level(const Range &ents, const BitRefLevel &bit, const bool only_tets=true, int verb=-1)
seed entities in the range and their adjacencies in a particular BitRefLevel
Definition: DeprecatedCoreInterface.cpp:37
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:95
MoFEM::Mat_Elastic_TransIso
Transverse Isotropic material data structure.
Definition: MaterialBlocks.hpp:384
MoFEM::Mat_Elastic_TransIso::data
_data_ data
Definition: MaterialBlocks.hpp:397
MoFEM::Mat_Interf
Linear interface data structure.
Definition: MaterialBlocks.hpp:443
MoFEM::Mat_Interf::data
_data_ data
Definition: MaterialBlocks.hpp:456
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:36
MoFEM::SnesCtx
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:27
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx
Definition: NonLinearElasticElement.hpp:570
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:40
NonlinearElasticElement::feRhs
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
Definition: NonLinearElasticElement.hpp:78
NonlinearElasticElement::setOfBlocks
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
Definition: NonLinearElasticElement.hpp:113
PostProcTemplateOnRefineMesh::postProcMesh
moab::Interface & postProcMesh
Definition: PostProcOnRefMesh.hpp:134
PostProcTemplateVolumeOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Definition: PostProcOnRefMesh.hpp:312
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:968
SmallStrainParaboloidalPlasticity
Small strain plasticity with paraboloidal yield criterion (Isotropic Hardening)
Definition: SmallStrainPlasticityMaterialModels.hpp:162
SmallStrainPlasticity
Small strain finite element implementation.
Definition: SmallStrainPlasticity.hpp:320
SmallStrainPlasticity::feRhs
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
Definition: SmallStrainPlasticity.hpp:360
SmallStrainPlasticity::createInternalVariablesTag
PetscErrorCode createInternalVariablesTag()
Definition: SmallStrainPlasticity.cpp:34
TimeForceScale
Force scale operator for reading two columns.
Definition: TimeForceScale.hpp:30
TimeForceScale::getForceScale
MoFEMErrorCode getForceScale(const double ts_t, double &scale)
Definition: TimeForceScale.hpp:101
Tri_Hand_tag
Definition: homogenisation_periodic_atom_test.cpp:71
VolumeCalculation
Calculate volume.
Definition: VolumeCalculation.hpp:31
xcoord_tag
Definition: homogenisation_periodic_atom_test.cpp:68
ycoord_tag
Definition: homogenisation_periodic_atom_test.cpp:69
zcoord_tag
Definition: homogenisation_periodic_atom_test.cpp:70
Generated on Tue Jul 5 2022 22:52:12 for MoFEM by Doxygen 1.9.4 and hosted at University of Glasgow