v0.13.1
rve_mech_plas_interface_disp.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 <BCs_RVELagrange_Disp.hpp>
39#include <VolumeCalculation.hpp>
40
41#include <CohesiveInterfaceElement.hpp>
42using namespace CohesiveElement;
43
44#ifndef WITH_ADOL_C
45#error "MoFEM need to be compiled with ADOL-C"
46#endif
47#include <adolc/adolc.h>
48#include <NonLinearElasticElement.hpp>
49#include <SmallTransverselyIsotropic.hpp>
50
51#include <adolc/adolc.h>
52#include <SmallStrainPlasticity.hpp>
53#include <SmallStrainPlasticityMaterialModels.hpp>
54
55
56#include <Projection10NodeCoordsOnField.hpp>
57#include <PostProcOnRefMesh.hpp>
58
59#include <string>
60
61using namespace boost::numeric;
62
63
64
65
66static char help[] =
67"...\n"
68"\n";
69
70int main(int argc, char *argv[]) {
71 MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
72
73 try {
74
75 moab::Core mb_instance;
76 moab::Interface& moab = mb_instance;
77 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
78 if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
79
80 PetscBool flg = PETSC_TRUE;
81 char mesh_file_name[255];
82 ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
83 if(flg != PETSC_TRUE) {
84 SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
85 }
86
87 PetscInt order;
88 ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
89 if(flg != PETSC_TRUE) {
90 order = 2;
91 }
92 PetscInt bubble_order;
93 ierr = PetscOptionsGetInt(PETSC_NULL,"-my_bubble_order",&bubble_order,&flg); CHKERRQ(ierr);
94 if(flg != PETSC_TRUE) {
95 bubble_order = order;
96 }
97
98 // use this if your mesh is partitioned and you run code on parts,
99 // you can solve very big problems
100 PetscBool is_partitioned = PETSC_FALSE;
101 ierr = PetscOptionsGetBool(PETSC_NULL,"-my_is_partitioned",&is_partitioned,&flg); CHKERRQ(ierr);
102
103 //Applied strain on the RVE (vector of length 6) strain=[xx, yy, zz, xy, xz, zy]^T
104 double mygiven_strain[6];
105 int nmax=6;
106 ierr = PetscOptionsGetRealArray(PETSC_NULL,"-my_given_strain",mygiven_strain,&nmax,&flg); CHKERRQ(ierr);
107 VectorDouble given_strain;
108 given_strain.resize(6);
109 cblas_dcopy(6, &mygiven_strain[0], 1, &given_strain(0), 1);
110 cout<<"given_strain ="<<given_strain<<endl;
111
112 //Read mesh to MOAB
113 if(is_partitioned == PETSC_TRUE) {
114 const char *option;
115 option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION=PARALLEL_PARTITION;";
116 rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
117 } else {
118 const char *option;
119 option = "";
120 rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
121 }
122
123
124// double young_modulus = 1;
125// SmallStrainParaboloidalPlasticity cp;
126// {
127// cp.tAgs.resize(3);
128// cp.tAgs[0] = 3;
129// cp.tAgs[1] = 4;
130// cp.tAgs[2] = 5;
131// cp.tOl = 1e-12;
132//
133// double poisson_ratio = 0.25;
134// cp.sIgma_yt = 1;
135// cp.sIgma_yc = 1;
136//
137// cp.Ht = 0.1;
138// cp.Hc = 0.1;
139//
140// cp.nup = 0.3;
141//
142//
143// {
144// ierr = PetscOptionsGetReal(0,"-my_young_modulus",&young_modulus,0); CHKERRQ(ierr);
145// ierr = PetscOptionsGetReal(0,"-my_poisson_ratio",&poisson_ratio,0); CHKERRQ(ierr);
146// cp.mu = MU(young_modulus,poisson_ratio);
147// cp.lambda = LAMBDA(young_modulus,poisson_ratio);
148// ierr = PetscOptionsGetReal(0,"-my_sigma_yt",&cp.sIgma_yt,0); CHKERRQ(ierr);
149// ierr = PetscOptionsGetReal(0,"-my_sigma_yc",&cp.sIgma_yc,0); CHKERRQ(ierr);
150//
151// ierr = PetscOptionsGetReal(0,"-my_Ht",&cp.Ht,0); CHKERRQ(ierr);
152// ierr = PetscOptionsGetReal(0,"-my_Hc",&cp.Hc,0); CHKERRQ(ierr);
153//
154// ierr = PetscOptionsGetReal(0,"-my_nup",&cp.nup,0); CHKERRQ(ierr);
155// }
156//
157// PetscPrintf(PETSC_COMM_WORLD,"young_modulus = %4.2e \n",young_modulus);
158// PetscPrintf(PETSC_COMM_WORLD,"poisson_ratio = %4.2e \n",poisson_ratio);
159//
160// PetscPrintf(PETSC_COMM_WORLD,"sIgma_yt = %4.2e \n",cp.sIgma_yt);
161// PetscPrintf(PETSC_COMM_WORLD,"sIgma_yc = %4.2e \n",cp.sIgma_yt);
162//
163// PetscPrintf(PETSC_COMM_WORLD,"nup = %4.2e \n",cp.nup);
164//
165// cp.sTrain.resize(6,false);
166// cp.sTrain.clear();
167// cp.plasticStrain.resize(6,false);
168// cp.plasticStrain.clear();
169// cp.internalVariables.resize(2,false);
170// cp.internalVariables.clear();
171// cp.createMatAVecR();
172// cp.snesCreate();
173// // ierr = SNESSetFromOptions(cp.sNes); CHKERRQ(ierr);
174// }
175
176
177 double young_modulus = 1;
178 SmallStrainJ2Plasticity cp;
179 {
180 cp.tAgs.resize(3);
181 cp.tAgs[0] = 3;
182 cp.tAgs[1] = 4;
183 cp.tAgs[2] = 5;
184 cp.tOl = 1e-12;
185
186 double poisson_ratio = 0.25;
187 cp.sIgma_y = 1;
188 cp.H = 0.1;
189 cp.K = 0;
190 cp.phi = 1;
191 {
192 ierr = PetscOptionsGetReal(0,"-my_young_modulus",&young_modulus,0); CHKERRQ(ierr);
193 ierr = PetscOptionsGetReal(0,"-my_poisson_ratio",&poisson_ratio,0); CHKERRQ(ierr);
194 cp.mu = MU(young_modulus,poisson_ratio);
195 cp.lambda = LAMBDA(young_modulus,poisson_ratio);
196 ierr = PetscOptionsGetReal(0,"-my_sigma_y",&cp.sIgma_y,0); CHKERRQ(ierr);
197 ierr = PetscOptionsGetReal(0,"-my_H",&cp.H,0); CHKERRQ(ierr);
198 ierr = PetscOptionsGetReal(0,"-my_K",&cp.K,0); CHKERRQ(ierr);
199 ierr = PetscOptionsGetReal(0,"-my_phi",&cp.phi,0); CHKERRQ(ierr);
200 }
201
202 PetscPrintf(PETSC_COMM_WORLD,"young_modulus = %4.2e \n",young_modulus);
203 PetscPrintf(PETSC_COMM_WORLD,"poisson_ratio = %4.2e \n",poisson_ratio);
204 PetscPrintf(PETSC_COMM_WORLD,"sIgma_y = %4.2e \n",cp.sIgma_y);
205 PetscPrintf(PETSC_COMM_WORLD,"H = %4.2e \n",cp.H);
206 PetscPrintf(PETSC_COMM_WORLD,"K = %4.2e \n",cp.K);
207 PetscPrintf(PETSC_COMM_WORLD,"phi = %4.2e \n",cp.phi);
208
209
210 cp.sTrain.resize(6,false);
211 cp.sTrain.clear();
212 cp.plasticStrain.resize(6,false);
213 cp.plasticStrain.clear();
214 cp.internalVariables.resize(7,false);
215 cp.internalVariables.clear();
216 cp.createMatAVecR();
217 cp.snesCreate();
218 // ierr = SNESSetFromOptions(cp.sNes); CHKERRQ(ierr);
219 }
220
221
222
223
224
225
226 //Create MoFEM (Joseph) database
227 MoFEM::Core core(moab);
228 MoFEM::Interface& m_field = core;
229
230 vector<BitRefLevel> bit_levels;
231 {
232 Tag th_meshset_info;
233 int def_meshset_info[2] = {0,0};
234 rval = moab.tag_get_handle(
235 "MESHSET_INFO",2,MB_TYPE_INTEGER,th_meshset_info,MB_TAG_CREAT|MB_TAG_SPARSE,&def_meshset_info
236 );
237 int meshset_data[2];
238 EntityHandle root = moab.get_root_set();
239 rval = moab.tag_get_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
240 if(meshset_data[0]==0) {
241 meshset_data[0] = 1;
242 rval = moab.tag_set_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
243
244 }
245 bit_levels.push_back(BitRefLevel().set(meshset_data[0]-1));
246 }
247
248
249 EntityHandle out_meshset;
250 rval = moab.create_meshset(MESHSET_SET,out_meshset); CHKERRQ_MOAB(rval);
251 // ierr = m_field.get_problem_finite_elements_entities("POTENTIAL_PROBLEM","POTENTIAL_ELEM",out_meshset); CHKERRQ(ierr);
252 ierr = m_field.get_entities_by_ref_level(bit_levels.back(),BitRefLevel().set(),out_meshset); CHKERRQ(ierr);
253 Range LatestRefinedTets;
254 rval = moab.get_entities_by_type(out_meshset, MBTET,LatestRefinedTets,true); CHKERRQ_MOAB(rval);
255 Range LatestRefinedPrisms;
256 rval = moab.get_entities_by_type(out_meshset, MBPRISM,LatestRefinedPrisms,true); CHKERRQ_MOAB(rval);
257
258 cout<< "========================== LatestRefinedPrisms "<<LatestRefinedPrisms.size()<<endl;
259 cout<< "========================== LatestRefinedTets "<<LatestRefinedTets.size()<<endl;
260
261
262 BitRefLevel problem_bit_level = bit_levels.back();
263 ierr = m_field.seed_ref_level_3D(0,problem_bit_level); CHKERRQ(ierr);
264
265 // const clock_t begin_time = clock();
266 ierr = m_field.build_fields(); CHKERRQ(ierr);
267 ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
268
269
270
271 Range prims_on_problem_bit_level;
272 ierr = m_field.get_entities_by_type_and_ref_level(
273 problem_bit_level,BitRefLevel().set(),MBPRISM,prims_on_problem_bit_level
274 ); CHKERRQ(ierr);
275 Range tets_on_problem_bit_level;
276 ierr = m_field.get_entities_by_type_and_ref_level(
277 problem_bit_level,BitRefLevel().set(),MBTET,tets_on_problem_bit_level
278 ); CHKERRQ(ierr);
279
280 //to create meshset from range
281 EntityHandle meshset_prims_on_problem_bit_level;
282 rval = moab.create_meshset(MESHSET_SET,meshset_prims_on_problem_bit_level); CHKERRQ_MOAB(rval);
283 rval = moab.add_entities(meshset_prims_on_problem_bit_level,prims_on_problem_bit_level); CHKERRQ_MOAB(rval);
284 ierr = m_field.seed_ref_level_MESHSET(meshset_prims_on_problem_bit_level,BitRefLevel().set()); CHKERRQ(ierr);
285
286
287 //Fields
288 ierr = m_field.add_field("DISPLACEMENT",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
289 ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
290 ierr = m_field.add_field("LAGRANGE_MUL_DISP",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
291
292 //add entitities (by tets) to the field
293 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"DISPLACEMENT"); CHKERRQ(ierr);
294 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
295
296 for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
297 if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
298 Range tris;
299 rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
300 ierr = m_field.add_ents_to_field_by_type(tris,MBTRI,"LAGRANGE_MUL_DISP"); CHKERRQ(ierr);
301 }
302 }
303
304
305 //set app. order
306 ierr = m_field.set_field_order(0,MBTET,"DISPLACEMENT",bubble_order); CHKERRQ(ierr);
307 ierr = m_field.set_field_order(0,MBTRI,"DISPLACEMENT",order); CHKERRQ(ierr);
308 ierr = m_field.set_field_order(0,MBEDGE,"DISPLACEMENT",order); CHKERRQ(ierr);
309 ierr = m_field.set_field_order(0,MBVERTEX,"DISPLACEMENT",1); CHKERRQ(ierr);
310
311 ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
312 ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
313 ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
314 ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
315
316 ierr = m_field.set_field_order(0,MBTRI,"LAGRANGE_MUL_DISP",order); CHKERRQ(ierr);
317 ierr = m_field.set_field_order(0,MBEDGE,"LAGRANGE_MUL_DISP",order); CHKERRQ(ierr);
318 ierr = m_field.set_field_order(0,MBVERTEX,"LAGRANGE_MUL_DISP",1); CHKERRQ(ierr);
319
320 //build field
321 ierr = m_field.build_fields(); CHKERRQ(ierr);
322
323 //10 node tets
324 Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
325 ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material,0); CHKERRQ(ierr);
326
327 //FE
328 ierr = m_field.add_finite_element("PLASTIC"); CHKERRQ(ierr);
329 //Define rows/cols and element data
330 ierr = m_field.modify_finite_element_add_field_row("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
331 ierr = m_field.modify_finite_element_add_field_col("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
332 ierr = m_field.modify_finite_element_add_field_data("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
333 ierr = m_field.modify_finite_element_add_field_data("PLASTIC","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
334
335
336
337 Range newtets;
338 for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,BLOCKSET|UNKNOWNNAME,it)) {
339 if(it->getName() != "MAT_PLASTIC") continue;
340 EntityHandle meshset = it->getMeshset();
341 int id = it->getMeshsetId();
342
343 rval = m_field.get_moab().get_entities_by_type(
344 meshset,MBTET,newtets,true
345 ); CHKERRQ_MOAB(rval);
346
347 //intersection of new and old tets for plastic
348 newtets = intersect(newtets,LatestRefinedTets);
349 ierr = m_field.seed_finite_elements(newtets); CHKERRQ(ierr);
350 }
351
352 ierr = m_field.add_ents_to_finite_element_by_type(newtets,MBTET,"PLASTIC"); CHKERRQ(ierr);
353// cout<< "========================== newtets "<<newtets.size()<<endl;
354
355
356 //================================================================================================================================
357 // Trans-Isotropic Yarns
358 //================================================================================================================================
359
360 NonlinearElasticElement trans_elastic(m_field,2);
361 trans_elastic.commonData.spatialPositions = "DISPLACEMENT";
362 trans_elastic.commonData.meshPositions = "MESH_NODE_POSITIONS";
363 std::map<int,boost::shared_ptr<SmallStrainTranverslyIsotropicADouble> > tranversly_isotropic_adouble_ptr_map;
364 std::map<int,boost::shared_ptr<SmallStrainTranverslyIsotropicDouble> > tranversly_isotropic_double_ptr_map;
365 bool trans_iso_blocks = false;
366 for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
367 //Get block name
368 string name = it->getName();
369 if (name.compare(0,20,"MAT_ELASTIC_TRANSISO") == 0) {
370 cout<<"================================ it is MAT_ELASTIC_TRANSISO "<<endl;
371 trans_iso_blocks = true;
372 int id = it->getMeshsetId();
373 Mat_Elastic_TransIso mydata;
374 ierr = it->getAttributeDataStructure(mydata); CHKERRQ(ierr);
375 tranversly_isotropic_adouble_ptr_map[id] = boost::make_shared<SmallStrainTranverslyIsotropicADouble>();
376 tranversly_isotropic_double_ptr_map[id] = boost::make_shared<SmallStrainTranverslyIsotropicDouble>();
377 //nu_p, nu_pz, E_p, E_z, G_zp
378 tranversly_isotropic_adouble_ptr_map.at(id)->E_p = mydata.data.Youngp;
379 tranversly_isotropic_double_ptr_map.at(id)->E_p = mydata.data.Youngp;
380 tranversly_isotropic_adouble_ptr_map.at(id)->E_z = mydata.data.Youngz;
381 tranversly_isotropic_double_ptr_map.at(id)->E_z = mydata.data.Youngz;
382 tranversly_isotropic_adouble_ptr_map.at(id)->nu_p = mydata.data.Poissonp;
383 tranversly_isotropic_double_ptr_map.at(id)->nu_p = mydata.data.Poissonp;
384 tranversly_isotropic_adouble_ptr_map.at(id)->nu_pz = mydata.data.Poissonpz;
385 tranversly_isotropic_double_ptr_map.at(id)->nu_pz = mydata.data.Poissonpz;
386 double shear_zp;
387 if(mydata.data.Shearzp!=0) {
388 shear_zp = mydata.data.Shearzp;
389 } else {
390 shear_zp = mydata.data.Youngz/(2*(1+mydata.data.Poissonpz));
391 }
392 tranversly_isotropic_adouble_ptr_map.at(it->getMeshsetId())->G_zp = shear_zp;
393 tranversly_isotropic_double_ptr_map.at(it->getMeshsetId())->G_zp = shear_zp;
394 //get tets from block where material is defined
395 EntityHandle meshset = it->getMeshset();
396 rval = m_field.get_moab().get_entities_by_type(
397 meshset,MBTET,trans_elastic.setOfBlocks[id].tEts,true
398 ); CHKERRQ_MOAB(rval);
399
400 //cout<<"================================ sit->second.tEts "<<sit->second.tEts.size()<<endl;
401 //intersection of new and old tets for tran-iso part
402 trans_elastic.setOfBlocks[id].tEts = intersect(trans_elastic.setOfBlocks[id].tEts,LatestRefinedTets);
403 //cout<<"================================ sit->second.tEts "<<sit->second.tEts.size()<<endl;
404
405 //adding material to nonlinear class
406 trans_elastic.setOfBlocks[id].iD = id;
407 //note that material parameters are defined internally in material model
408 trans_elastic.setOfBlocks[id].E = 0; // this is not working for this material
409 trans_elastic.setOfBlocks[id].PoissonRatio = 0; // this is not working for this material
410 trans_elastic.setOfBlocks[id].materialDoublePtr = tranversly_isotropic_double_ptr_map.at(id);
411 trans_elastic.setOfBlocks[id].materialAdoublePtr = tranversly_isotropic_adouble_ptr_map.at(id);
412 ierr = m_field.seed_finite_elements(trans_elastic.setOfBlocks[id].tEts); CHKERRQ(ierr);
413 }
414 }
415
416
417
418 ierr = m_field.add_finite_element("TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
419 ierr = m_field.add_finite_element("TRAN_ISOTROPIC_ELASTIC",MF_ZERO); CHKERRQ(ierr);
420 ierr = m_field.modify_finite_element_add_field_row("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
421 ierr = m_field.modify_finite_element_add_field_col("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
422 ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
423 ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","POTENTIAL_FIELD"); CHKERRQ(ierr);
424 if(m_field.check_field("MESH_NODE_POSITIONS")) {
425 ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
426 }
427 for(
428 map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
429 sit!=trans_elastic.setOfBlocks.end();sit++
430 ) {
431 ierr = m_field.add_ents_to_finite_element_by_type(sit->second.tEts,MBTET,"TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
432 }
433
434
435 //Rhs
436 trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("DISPLACEMENT",trans_elastic.commonData));
437 trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("POTENTIAL_FIELD",trans_elastic.commonData));
438 if(m_field.check_field("MESH_NODE_POSITIONS")) {
439 trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("MESH_NODE_POSITIONS",trans_elastic.commonData));
440 }
441 map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
442 for(;sit!=trans_elastic.setOfBlocks.end();sit++) {
443 trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress("DISPLACEMENT",sit->second,trans_elastic.commonData,2,false,false,true));
444 trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpRhsPiolaKirchhoff("DISPLACEMENT",sit->second,trans_elastic.commonData));
445 }
446
447 //Lhs
448 trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("DISPLACEMENT",trans_elastic.commonData));
449 trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("POTENTIAL_FIELD",trans_elastic.commonData));
450 if(m_field.check_field("MESH_NODE_POSITIONS")) {
451 trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("MESH_NODE_POSITIONS",trans_elastic.commonData));
452 }
453 sit = trans_elastic.setOfBlocks.begin();
454 for(;sit!=trans_elastic.setOfBlocks.end();sit++) {
455 trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress("DISPLACEMENT",sit->second,trans_elastic.commonData,2,true,false,true));
456 trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpLhsPiolaKirchhoff_dx("DISPLACEMENT","DISPLACEMENT",sit->second,trans_elastic.commonData));
457 }
458
459 //================================================================================================================================
460 // Interface Cohesive elemetns
461 //================================================================================================================================
462 //FE Interface
463 ierr = m_field.add_finite_element("INTERFACE"); CHKERRQ(ierr);
464 ierr = m_field.modify_finite_element_add_field_row("INTERFACE","DISPLACEMENT"); CHKERRQ(ierr);
465 ierr = m_field.modify_finite_element_add_field_col("INTERFACE","DISPLACEMENT"); CHKERRQ(ierr);
466 ierr = m_field.modify_finite_element_add_field_data("INTERFACE","DISPLACEMENT"); CHKERRQ(ierr);
467 ierr = m_field.modify_finite_element_add_field_data("INTERFACE","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
468
469 ierr = m_field.add_ents_to_finite_element_by_bit_ref(problem_bit_level,BitRefLevel().set(),"INTERFACE",MBPRISM); CHKERRQ(ierr);
470
471 boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation> interface_materials;
472
473 //FIXME this in fact allow only for one type of interface,
474 //problem is Young Moduls in interface mayoung_modulusterial
475
476
477 double interface_beta, interface_ft, interface_Gf, interface_h;
478 ierr = PetscOptionsGetReal(0,"-interface_beta",&interface_beta,0); CHKERRQ(ierr);
479 ierr = PetscOptionsGetReal(0,"-interface_ft",&interface_ft,0); CHKERRQ(ierr);
480 ierr = PetscOptionsGetReal(0,"-interface_Gf",&interface_Gf,0); CHKERRQ(ierr);
481 ierr = PetscOptionsGetReal(0,"-interface_h",&interface_h,0); CHKERRQ(ierr);
482
483 for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
484 cout << endl << *it << endl;
485 //Get block name
486 string name = it->getName();
487 if (name.compare(0,10,"MAT_INTERF") == 0) {
488 Mat_Interf mydata;
489 ierr = it->getAttributeDataStructure(mydata); CHKERRQ(ierr);
490 cout << mydata;
491 interface_materials.push_back(new CohesiveInterfaceElement::PhysicalEquation(m_field));
492// interface_materials.back().h = 0.134; // we used this because E0=35Gpa and Em=4.7Gpa
493 interface_materials.back().youngModulus = mydata.data.alpha;
494// interface_materials.back().beta = mydata.data.beta;
495// interface_materials.back().ft = mydata.data.ft;
496// interface_materials.back().Gf = mydata.data.Gf;
497
498 interface_materials.back().h = interface_h;
499 interface_materials.back().beta = interface_beta;
500 interface_materials.back().ft = interface_ft;
501 interface_materials.back().Gf = interface_Gf;
502
503 EntityHandle meshset = it->getMeshset();
504 Range tris;
505 rval = moab.get_entities_by_type(meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
506 Range ents3d;
507 rval = moab.get_adjacencies(tris,3,false,ents3d,moab::Interface::UNION); CHKERRQ_MOAB(rval);
508 interface_materials.back().pRisms = ents3d.subset_by_type(MBPRISM);
509 }
510 }
511
512 // cout<<"young_modulus "<<young_modulus<<endl;
513 { //FIXME
514 boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>::iterator pit = interface_materials.begin();
515 for(; pit != interface_materials.end();pit++) {
516 pit->youngModulus = young_modulus; //Young's modulus of bulk elastic material (then for interface E0=E/h)
517 }
518 }
519 CohesiveInterfaceElement cohesive_elements(m_field);
520 ierr = cohesive_elements.addOps("DISPLACEMENT",interface_materials); CHKERRQ(ierr);
521 //================================================================================================================================
522
523 BCs_RVELagrange_Disp lagrangian_element_disp(m_field);
524 lagrangian_element_disp.addLagrangiangElement("LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS");
525
526
527
528
529 //build finite elements
530 ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
531 //build adjacencies
532 ierr = m_field.build_adjacencies(problem_bit_level); CHKERRQ(ierr);
533
534
535 DMType dm_name = "PLASTIC_PROB";
536 ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
537
538 DM dm;
539 ierr = DMCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(ierr);
540 ierr = DMSetType(dm,dm_name);CHKERRQ(ierr);
541
542 //set dm datastruture which created mofem datastructures
543 ierr = DMMoFEMCreateMoFEM(dm,&m_field,dm_name,problem_bit_level); CHKERRQ(ierr);
544 ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
545 ierr = DMMoFEMSetIsPartitioned(dm,is_partitioned); CHKERRQ(ierr);
546 //add elements to dm
547 ierr = DMMoFEMAddElement(dm,"PLASTIC"); CHKERRQ(ierr);
548 ierr = DMMoFEMAddElement(dm,"LAGRANGE_ELEM"); CHKERRQ(ierr);
549 ierr = DMMoFEMAddElement(dm,"TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
550 ierr = DMMoFEMAddElement(dm,"INTERFACE"); CHKERRQ(ierr);
551
552 ierr = DMSetUp(dm); CHKERRQ(ierr);
553
554 //create matrices
555 Vec F,D;
556 Mat Aij;
557 ierr = DMCreateGlobalVector_MoFEM(dm,&D); CHKERRQ(ierr);
558 ierr = VecDuplicate(D,&F); CHKERRQ(ierr);
559 ierr = DMCreateMatrix_MoFEM(dm,&Aij); CHKERRQ(ierr);
560 ierr = VecZeroEntries(D); CHKERRQ(ierr);
561 ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
562 ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
563 ierr = DMoFEMMeshToLocalVector(dm,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
564 ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
565
566 vector<Vec> Fvec(6); //jthis will be used to caluclate the homogenised stiffness matix
567 for(int ii = 0;ii<6;ii++) {
568 ierr = VecDuplicate(D,&Fvec[ii]); CHKERRQ(ierr);
569 ierr = VecZeroEntries(Fvec[ii]); CHKERRQ(ierr);
570 ierr = VecGhostUpdateBegin(Fvec[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
571 ierr = VecGhostUpdateEnd(Fvec[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
572 }
573
574 //assemble Aij and F
575 SmallStrainPlasticity small_strain_plasticity(m_field);
576 {
577 PetscBool bbar = PETSC_TRUE;
578 ierr = PetscOptionsGetBool(0,"-my_bbar",&bbar,0); CHKERRQ(ierr);
579 small_strain_plasticity.commonData.bBar = bbar;
580 }
581 {
582 small_strain_plasticity.feRhs.getOpPtrVector().push_back(
583 new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
584 "DISPLACEMENT",small_strain_plasticity.commonData
585 )
586 );
587 small_strain_plasticity.feRhs.getOpPtrVector().push_back(
588 new SmallStrainPlasticity::OpCalculateStress(
589 m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp
590 )
591 );
592 small_strain_plasticity.feRhs.getOpPtrVector().push_back(
593 new SmallStrainPlasticity::OpAssembleRhs(
594 "DISPLACEMENT",small_strain_plasticity.commonData
595 )
596 );
597 small_strain_plasticity.feLhs.getOpPtrVector().push_back(
598 new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
599 "DISPLACEMENT",small_strain_plasticity.commonData
600 )
601 );
602 small_strain_plasticity.feLhs.getOpPtrVector().push_back(
603 new SmallStrainPlasticity::OpCalculateStress(
604 m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp
605 )
606 );
607 small_strain_plasticity.feLhs.getOpPtrVector().push_back(
608 new SmallStrainPlasticity::OpAssembleLhs(
609 "DISPLACEMENT",small_strain_plasticity.commonData
610 )
611 );
612 small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
613 new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
614 "DISPLACEMENT",small_strain_plasticity.commonData
615 )
616 );
617 small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
618 new SmallStrainPlasticity::OpCalculateStress(
619 m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp
620 )
621 );
622 small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
623 new SmallStrainPlasticity::OpUpdate(
624 "DISPLACEMENT",small_strain_plasticity.commonData
625 )
626 );
627 ierr = small_strain_plasticity.createInternalVariablesTag(); CHKERRQ(ierr);
628 }
629
630 lagrangian_element_disp.setRVEBCsOperatorsNonlinear("DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS",Aij,Fvec,F,given_strain);
631 TimeForceScale time_force_scale("-my_load_history",false);
632 lagrangian_element_disp.methodsOp.push_back(new TimeForceScale("-my_macro_strian_history",false));
633
634
635
636 //Adding elements to DMSnes
637 //Rhs
638 ierr = DMMoFEMSNESSetFunction(dm,"PLASTIC",&small_strain_plasticity.feRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
639 ierr = DMMoFEMSNESSetFunction(dm,"TRAN_ISOTROPIC_ELASTIC",&trans_elastic.feRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
640 ierr = DMMoFEMSNESSetFunction(dm,"INTERFACE",&cohesive_elements.feRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
641 ierr = DMMoFEMSNESSetFunction(dm,"LAGRANGE_ELEM",&lagrangian_element_disp.feRVEBCRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
642 ierr = DMMoFEMSNESSetFunction(dm,"LAGRANGE_ELEM",&lagrangian_element_disp.feRVEBCRhsResidual,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
643
644 //Lhs
645 ierr = DMMoFEMSNESSetJacobian(dm,"PLASTIC",&small_strain_plasticity.feLhs,NULL,NULL); CHKERRQ(ierr);
646 ierr = DMMoFEMSNESSetJacobian(dm,"TRAN_ISOTROPIC_ELASTIC",&trans_elastic.feLhs,NULL,NULL); CHKERRQ(ierr);
647 ierr = DMMoFEMSNESSetJacobian(dm,"INTERFACE",&cohesive_elements.feLhs,NULL,NULL); CHKERRQ(ierr);
648 ierr = DMMoFEMSNESSetJacobian(dm,"LAGRANGE_ELEM",&lagrangian_element_disp.feRVEBCLhs,NULL,NULL); CHKERRQ(ierr);
649
650 // Create Time Stepping solver
651 SNES snes;
652 SnesCtx *snes_ctx;
653
654 ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
655 //ierr = SNESSetDM(snes,dm); CHKERRQ(ierr);
656 ierr = DMMoFEMGetSnesCtx(dm,&snes_ctx); CHKERRQ(ierr);
657 ierr = SNESSetFunction(snes,F,SnesRhs,snes_ctx); CHKERRQ(ierr);
658 ierr = SNESSetJacobian(snes,Aij,Aij,SnesMat,snes_ctx); CHKERRQ(ierr);
659 ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
660
661
662 PostProcVolumeOnRefinedMesh post_proc(m_field);
663 ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
664 ierr = post_proc.addFieldValuesPostProc("DISPLACEMENT"); CHKERRQ(ierr);
665 ierr = post_proc.addFieldValuesGradientPostProc("DISPLACEMENT"); CHKERRQ(ierr);
666 ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
667
668 // Volume calculation
669 //==============================================================================================================================
670 VolumeElementForcesAndSourcesCore vol_elem(m_field);
671 Vec volume_vec;
672 int volume_vec_ghost[] = { 0 };
673 ierr = VecCreateGhost(
674 PETSC_COMM_WORLD,(!m_field.get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
675 ); CHKERRQ(ierr);
676 ierr = VecZeroEntries(volume_vec); CHKERRQ(ierr);
677 vol_elem.getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
678
679 ierr = m_field.loop_finite_elements("PLASTIC_PROB","PLASTIC", vol_elem); CHKERRQ(ierr);
680 ierr = m_field.loop_finite_elements("PLASTIC_PROB","TRAN_ISOTROPIC_ELASTIC", vol_elem); CHKERRQ(ierr);
681
682 ierr = VecAssemblyBegin(volume_vec); CHKERRQ(ierr);
683 ierr = VecAssemblyEnd(volume_vec); CHKERRQ(ierr);
684 double rve_volume;
685 ierr = VecSum(volume_vec,&rve_volume); CHKERRQ(ierr);
686
687 // ierr = VecView(volume_vec,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
688 ierr = PetscPrintf(PETSC_COMM_WORLD,"RVE Volume %3.2g\n",rve_volume); CHKERRQ(ierr);
689 ierr = VecDestroy(&volume_vec);
690 //=============================================================================================================================
691
692 double final_time = 1,delta_time = 0.1;
693 ierr = PetscOptionsGetReal(0,"-my_final_time",&final_time,0); CHKERRQ(ierr);
694 ierr = PetscOptionsGetReal(0,"-my_delta_time",&delta_time,0); CHKERRQ(ierr);
695 double delta_time0 = delta_time;
696
697// ierr = MatView(Aij,PETSC_VIEWER_DRAW_SELF); CHKERRQ(ierr);
698// std::string wait;
699// std::cin >> wait;
700
701 Vec D0;
702 ierr = VecDuplicate(D,&D0); CHKERRQ(ierr);
703
704 int step = 0;
705 double t = 0;
706 SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
707 for(;t<final_time;step++) {
708 t += delta_time;
709
710 if(t>final_time){
711 break;
712 }
713 PetscPrintf(
714 PETSC_COMM_WORLD,"Step %d Time %6.4g final time %3.2g\n",step,t,final_time
715 );
716 //set time
717 lagrangian_element_disp.getLoopFeRVEBCRhs().ts_t = t;
718 //solve problem at step
719 ierr = VecAssemblyBegin(D);
720 ierr = VecAssemblyEnd(D);
721 ierr = VecCopy(D,D0); CHKERRQ(ierr);
722 if(step == 0 || reason < 0) {
723 ierr = SNESSetLagJacobian(snes,-2); CHKERRQ(ierr);
724 } else {
725 ierr = SNESSetLagJacobian(snes,-1); CHKERRQ(ierr);
726 }
727 ierr = SNESSolve(snes,PETSC_NULL,D); CHKERRQ(ierr);
728
729 int its;
730 ierr = SNESGetIterationNumber(snes,&its); CHKERRQ(ierr);
731
732 ierr = PetscPrintf(
733 PETSC_COMM_WORLD,"number of Newton iterations = %D\n",its
734 ); CHKERRQ(ierr);
735
736 ierr = SNESGetConvergedReason(snes,&reason); CHKERRQ(ierr);
737
738 if(reason<0) {
739 t -= delta_time;
740 delta_time *= 0.1;
741 ierr = VecCopy(D0,D); CHKERRQ(ierr);
742 } else {const int its_d = 60;
743 const double gamma = 0.5;
744 const double max_reudction = 1;
745 const double min_reduction = 1e-1;
746 double reduction;
747 reduction = pow((double)its_d/(double)(its+1),gamma);
748 if(delta_time >= max_reudction*delta_time0 && reduction > 1) {
749 reduction = 1;
750 } else if(delta_time <= min_reduction*delta_time0 && reduction < 1) {
751 reduction = 1;
752 }
753
754 ierr = PetscPrintf(
755 PETSC_COMM_WORLD,
756 "reduction delta_time = %6.4e delta_time = %6.4e\n",
757 reduction,delta_time
758 ); CHKERRQ(ierr);
759 delta_time *= reduction;
760 if(reduction>1 && delta_time < min_reduction*delta_time0) {
761 delta_time = min_reduction*delta_time0;
762 }
763
764 ierr = DMoFEMMeshToGlobalVector(
765 dm,D,INSERT_VALUES,SCATTER_REVERSE
766 ); CHKERRQ(ierr);
767
768 ierr = DMoFEMLoopFiniteElements(
769 dm,"PLASTIC",&small_strain_plasticity.feUpdate
770 ); CHKERRQ(ierr);
771
772
773
774 ierr = DMoFEMLoopFiniteElements(
775 dm,"INTERFACE",&cohesive_elements.feHistory
776 ); CHKERRQ(ierr);
777
778
779 //Save data on mesh
780 {
781 ierr = DMoFEMLoopFiniteElements(
782 dm,"PLASTIC",&post_proc
783 ); CHKERRQ(ierr);
784 string out_file_name;
785 {
786 std::ostringstream stm;
787 stm << "out_" << step << ".h5m";
788 out_file_name = stm.str();
789 }
790 ierr = PetscPrintf(
791 PETSC_COMM_WORLD,"out file %s\n",out_file_name.c_str()
792 ); CHKERRQ(ierr);
793 rval = post_proc.postProcMesh.write_file(
794 out_file_name.c_str(),"MOAB","PARALLEL=WRITE_PART"
795 ); CHKERRQ_MOAB(rval);
796 }
797
798 //===================================== Homgenised stress (for given strain) ====================================================
799 VectorDouble StressHomo;
800 StressHomo.resize(6);
801 StressHomo.clear();
802
803 //calculate honmogenised stress
804 //create a vector for 6 components of homogenized stress
805 Vec stress_homo;
806 int stress_homo_ghost[] = { 0,1,2,3,4,5,6 };
807 ierr = VecCreateGhost(
808 PETSC_COMM_WORLD,(!m_field.get_comm_rank())?6:0,6,6,stress_homo_ghost,&stress_homo
809 ); CHKERRQ(ierr);
810
811 lagrangian_element_disp.setRVEBCsHomoStressOperators("DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS",stress_homo);
812
813
814 PetscScalar *avec;
815 ierr = VecGetArray(stress_homo,&avec); CHKERRQ(ierr);
816 ierr = VecZeroEntries(stress_homo); CHKERRQ(ierr);
817 ierr = m_field.loop_finite_elements(
818 "PLASTIC_PROB","LAGRANGE_ELEM",lagrangian_element_disp.getLoopFeRVEBCStress()
819 ); CHKERRQ(ierr);
820 ierr = PetscOptionsGetReal(
821 PETSC_NULL,"-my_rve_volume",&rve_volume,PETSC_NULL
822 ); CHKERRQ(ierr);
823 ierr = VecAssemblyBegin(stress_homo); CHKERRQ(ierr);
824 ierr = VecAssemblyEnd(stress_homo); CHKERRQ(ierr);
825 ierr = VecScale(stress_homo,1.0/rve_volume); CHKERRQ(ierr);
826 ierr = VecGhostUpdateBegin(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
827 ierr = VecGhostUpdateEnd(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
828
829 for(int jj=0; jj<6; jj++) {
830 StressHomo(jj) = avec[jj];
831 }
832
833 double scale;
834 ierr = time_force_scale.getForceScale(t,scale); CHKERRQ(ierr);
835
836 PetscPrintf(PETSC_COMM_WORLD,
837 "Macro_Strain Homo_Stress Path "
838 );
839
840 //cout<<"Macro-strain ";
841 for(int ii=0; ii<6; ii++) {
842 PetscPrintf(
843 PETSC_COMM_WORLD,
844 "%8.5e\t",
845 t*given_strain(ii)
846 );
847 }
848
849 //cout<<"Homo stress ";
850 for(int ii=0; ii<6; ii++) {
851 PetscPrintf(
852 PETSC_COMM_WORLD,
853 "%8.5e\t",
854 StressHomo(ii)
855 );
856 }
857
858 PetscPrintf(PETSC_COMM_WORLD,
859 "\n"
860 );
861// //==============================================================================================================================
862 }
863 }
864 ierr = VecDestroy(&D0); CHKERRQ(ierr);
865
866 ierr = MatDestroy(&Aij); CHKERRQ(ierr);
867 ierr = VecDestroy(&F); CHKERRQ(ierr);
868 ierr = VecDestroy(&D); CHKERRQ(ierr);
869 }
870 CATCH_ERRORS;
871
872 MoFEM::Core::Finalize();
873}
BCs_RVELagrange_Disp.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
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::add_ents_to_finite_element_by_bit_ref
virtual MoFEMErrorCode add_ents_to_finite_element_by_bit_ref(const BitRefLevel &bit, const BitRefLevel &mask, const std::string &name, EntityType type, int verb=DEFAULT_VERBOSITY)=0
add TET entities from given refinement level to finite element database given by name
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
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
main
int main(int argc, char *argv[])
Definition: rve_mech_plas_interface_disp.cpp:70
help
static char help[]
Definition: rve_mech_plas_interface_disp.cpp:66
BCs_RVELagrange_Disp
Definition: BCs_RVELagrange_Disp.hpp:18
BCs_RVELagrange_Disp::feRVEBCLhs
MyTriFE feRVEBCLhs
Definition: BCs_RVELagrange_Disp.hpp:30
BCs_RVELagrange_Disp::methodsOp
boost::ptr_vector< MethodForForceScaling > methodsOp
Definition: BCs_RVELagrange_Disp.hpp:25
BCs_RVELagrange_Disp::getLoopFeRVEBCRhs
MyTriFE & getLoopFeRVEBCRhs()
Definition: BCs_RVELagrange_Disp.hpp:36
BCs_RVELagrange_Disp::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_Disp.hpp:535
BCs_RVELagrange_Disp::feRVEBCRhsResidual
MyTriFE feRVEBCRhsResidual
Definition: BCs_RVELagrange_Disp.hpp:32
BCs_RVELagrange_Disp::setRVEBCsHomoStressOperators
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec Stress_Homo)
Definition: BCs_RVELagrange_Disp.hpp:671
BCs_RVELagrange_Disp::getLoopFeRVEBCStress
MyTriFE & getLoopFeRVEBCStress()
Definition: BCs_RVELagrange_Disp.hpp:39
BCs_RVELagrange_Disp::addLagrangiangElement
PetscErrorCode addLagrangiangElement(const string element_name, const string field_name, const string lagrang_field_name, const string mesh_nodals_positions)
Definition: BCs_RVELagrange_Disp.hpp:58
BCs_RVELagrange_Disp::feRVEBCRhs
MyTriFE feRVEBCRhs
Definition: BCs_RVELagrange_Disp.hpp:29
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
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_3D
DEPRECATED MoFEMErrorCode seed_ref_level_3D(const EntityHandle meshset, const BitRefLevel &bit, int verb=-1)
seed 2D entities in the meshset and their adjacencies (only TETs adjacencies) in a particular BitRefL...
Definition: DeprecatedCoreInterface.cpp:31
MoFEM::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
SmallStrainJ2Plasticity
J2 plasticity (Kinematic Isotropic (Linear) Hardening)
Definition: SmallStrainPlasticityMaterialModels.hpp:32
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
VolumeCalculation
Calculate volume.
Definition: VolumeCalculation.hpp:31
Generated on Tue Jul 5 2022 22:52:12 for MoFEM by Doxygen 1.9.4 and hosted at University of Glasgow