v0.13.2
Loading...
Searching...
No Matches
Public Member Functions | Public Attributes | List of all members
EshelbianPlasticity::EshelbianCore Struct Reference

#include <users_modules/eshelbian_plasticty/src/EshelbianPlasticity.hpp>

Inheritance diagram for EshelbianPlasticity::EshelbianCore:
[legend]
Collaboration diagram for EshelbianPlasticity::EshelbianCore:
[legend]

Public Member Functions

MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const
 Getting interface of core database. More...
 
 EshelbianCore (MoFEM::Interface &m_field)
 
virtual ~EshelbianCore ()
 
MoFEMErrorCode getOptions ()
 
template<typename BC >
MoFEMErrorCode getBc (boost::shared_ptr< BC > &bc_vec_ptr, const std::string block_set_name, const int nb_attributes)
 
MoFEMErrorCode getSpatialDispBc ()
 
MoFEMErrorCode getSpatialRotationBc ()
 
MoFEMErrorCode getSpatialTractionBc ()
 
MoFEMErrorCode getTractionFreeBc (const EntityHandle meshset, boost::shared_ptr< TractionFreeBc > &bc_ptr, const std::string disp_block_set_name, const std::string rot_block_set_name)
 
MoFEMErrorCode getSpatialTractionFreeBc (const EntityHandle meshset=0)
 
MoFEMErrorCode addFields (const EntityHandle meshset=0)
 
MoFEMErrorCode addVolumeFiniteElement (const EntityHandle meshset=0)
 
MoFEMErrorCode addBoundaryFiniteElement (const EntityHandle meshset=0)
 
MoFEMErrorCode addDMs (const BitRefLevel bit=BitRefLevel().set(0))
 
MoFEMErrorCode addMaterial_HMHHStVenantKirchhoff (const int tape, const double lambda, const double mu, const double sigma_y)
 
MoFEMErrorCode addMaterial_HMHMooneyRivlin (const int tape, const double alpha, const double beta, const double lambda, const double sigma_y)
 
MoFEMErrorCode setBaseVolumeElementOps (const int tag, const bool do_rhs, const bool do_lhs, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > &fe)
 
MoFEMErrorCode setGenericVolumeElementOps (const int tag, const bool add_elastic, const bool add_material, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > &fe_rhs, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > &fe_lhs)
 
MoFEMErrorCode setGenericFaceElementOps (const bool add_elastic, const bool add_material, boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > &fe_rhs, boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > &fe_lhs)
 
MoFEMErrorCode setElasticElementOps (const int tag)
 
MoFEMErrorCode setElasticElementToTs (DM dm)
 
MoFEMErrorCode setUpTSElastic (TS ts, Mat m, Vec f)
 
MoFEMErrorCode solveElastic (TS ts, Vec x)
 
MoFEMErrorCode postProcessResults (const int tag, const std::string file)
 
- Public Member Functions inherited from MoFEM::UnknownInterface
virtual MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const =0
 
template<class IFACE >
MoFEMErrorCode registerInterface (bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE *&iface) const
 Get interface refernce to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE **const iface) const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get reference to interface. More...
 
template<class IFACE >
IFACE * getInterface () const
 Function returning pointer to interface. More...
 
virtual ~UnknownInterface ()=default
 

Public Attributes

MoFEM::InterfacemField
 
boost::shared_ptr< DataAtIntegrationPtsdataAtPts
 
boost::shared_ptr< PhysicalEquationsphysicalEquations
 
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > elasticFeRhs
 
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > elasticFeLhs
 
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > elasticBcLhs
 
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > elasticBcRhs
 
boost::shared_ptr< EpFEMethodschurAssembly
 
SmartPetscObj< DM > dM
 Coupled problem all fields. More...
 
SmartPetscObj< DM > dmElastic
 Elastic problem. More...
 
SmartPetscObj< DM > dmElasticSchurStreach
 Sub problem of dmElastic Schur. More...
 
SmartPetscObj< DM > dmElasticSchurBubble
 Sub problem of dmElastic Schur. More...
 
SmartPetscObj< DM > dmElasticSchurSpatialDisp
 Sub problem of dmElastic Schur. More...
 
SmartPetscObj< DM > dmElasticSchurOmega
 Sub problem of dmElastic Schur. More...
 
SmartPetscObj< DM > dmMaterial
 Material problem. More...
 
const std::string piolaStress
 
const std::string eshelbyStress
 
const std::string spatialDisp
 
const std::string materialDisp
 
const std::string streachTensor
 
const std::string rotAxis
 
const std::string materialGradient
 
const std::string tauField
 
const std::string lambdaField
 
const std::string bubbleField
 
const std::string elementVolumeName
 
const std::string naturalBcElement
 
const std::string essentialBcElement
 
int spaceOrder
 
int materialOrder
 
double alphaU
 
double alphaW
 
double alphaRho
 
double precEps
 
boost::shared_ptr< BcDispVecbcSpatialDispVecPtr
 
boost::shared_ptr< BcRotVecbcSpatialRotationVecPtr
 
boost::shared_ptr< TractionBcVecbcSpatialTraction
 
boost::shared_ptr< TractionFreeBcbcSpatialFreeTraction
 

Additional Inherited Members

- Static Public Member Functions inherited from MoFEM::UnknownInterface
static MoFEMErrorCode getLibVersion (Version &version)
 Get library version. More...
 
static MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version)
 Get database major version. More...
 
static MoFEMErrorCode setFileVersion (moab::Interface &moab, Version version=Version(MoFEM_VERSION_MAJOR, MoFEM_VERSION_MINOR, MoFEM_VERSION_BUILD))
 Get database major version. More...
 
static MoFEMErrorCode getInterfaceVersion (Version &version)
 Get database major version. More...
 

Detailed Description

Examples
ep.cpp.

Definition at line 1157 of file EshelbianPlasticity.hpp.

Constructor & Destructor Documentation

◆ EshelbianCore()

EshelbianPlasticity::EshelbianCore::EshelbianCore ( MoFEM::Interface m_field)
Examples
EshelbianPlasticity.cpp.

Definition at line 40 of file EshelbianPlasticity.cpp.

41 : mField(m_field), piolaStress("P"), eshelbyStress("S"), spatialDisp("w"),
42 materialDisp("W"), streachTensor("u"), rotAxis("omega"),
43 materialGradient("G"), tauField("TAU"), lambdaField("LAMBDA"),
44 bubbleField("BUBBLE"), elementVolumeName("EP"),
45 naturalBcElement("NATURAL_BC"), essentialBcElement("ESSENTIAL_BC") {
46
47 ierr = getOptions();
48 CHKERRABORT(PETSC_COMM_WORLD, ierr);
49}
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76

◆ ~EshelbianCore()

EshelbianPlasticity::EshelbianCore::~EshelbianCore ( )
virtual
Examples
EshelbianPlasticity.cpp.

Definition at line 51 of file EshelbianPlasticity.cpp.

51{}

Member Function Documentation

◆ addBoundaryFiniteElement()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::addBoundaryFiniteElement ( const EntityHandle  meshset = 0)
Examples
EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 254 of file EshelbianPlasticity.cpp.

254 {
256
257 auto bc_elements_add_to_range = [&](const std::string disp_block_set_name,
258 Range &r) {
261 if (it->getName().compare(0, disp_block_set_name.length(),
262 disp_block_set_name) == 0) {
263 Range faces;
264 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
265 true);
266 r.merge(faces);
267 }
268 }
270 };
271
272 // set finite element fields
273 auto add_field_to_fe = [this](const std::string fe,
274 const std::string field_name) {
280 };
281
282 Range natural_bc_elements;
283 CHKERR bc_elements_add_to_range("SPATIAL_DISP_BC", natural_bc_elements);
284 CHKERR bc_elements_add_to_range("SPATIAL_ROTATION_BC", natural_bc_elements);
285 Range essentail_bc_elements;
286 CHKERR bc_elements_add_to_range("SPATIAL_TRACTION_BC", essentail_bc_elements);
287
289 CHKERR mField.add_ents_to_finite_element_by_type(natural_bc_elements, MBTRI,
291 CHKERR add_field_to_fe(naturalBcElement, piolaStress);
292 CHKERR add_field_to_fe(naturalBcElement, eshelbyStress);
294
296 CHKERR mField.add_ents_to_finite_element_by_type(essentail_bc_elements, MBTRI,
298 CHKERR add_field_to_fe(essentialBcElement, piolaStress);
299 CHKERR add_field_to_fe(essentialBcElement, eshelbyStress);
301
303}
@ MF_ZERO
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BLOCKSET
Definition: definitions.h:148
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
#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.
int r
Definition: sdf.py:5
constexpr auto field_name
virtual moab::Interface & get_moab()=0

◆ addDMs()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::addDMs ( const BitRefLevel  bit = BitRefLevel().set(0))
Examples
EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 305 of file EshelbianPlasticity.cpp.

305 {
307
308 // find adjacencies between finite elements and dofs
310
311 // Create coupled problem
312 dM = createSmartDM(mField.get_comm(), "DMMOFEM");
313 CHKERR DMMoFEMCreateMoFEM(dM, &mField, "ESHELBY_PLASTICITY", bit,
314 BitRefLevel().set());
316 CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
320 mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
321 CHKERR DMSetUp(dM);
322 mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
323
324 auto remove_dofs_on_essential_spatial_stress_boundary =
325 [&](const std::string prb_name) {
327 for (int d : {0, 1, 2})
329 prb_name, piolaStress, (*bcSpatialFreeTraction)[d], d, d, 0, 100,
330 NOISY, true);
332 };
333 CHKERR remove_dofs_on_essential_spatial_stress_boundary("ESHELBY_PLASTICITY");
334
335 // Create elastic sub-problem
336 dmElastic = createSmartDM(mField.get_comm(), "DMMOFEM");
337 CHKERR DMMoFEMCreateSubDM(dmElastic, dM, "ELASTIC_PROBLEM");
349 CHKERR DMSetUp(dmElastic);
350
351 // Create elastic streach-problem
354 "ELASTIC_PROBLEM_STREACH_SCHUR");
366
367 // Create elastic bubble-problem
370 "ELASTIC_PROBLEM_BUBBLE_SCHUR");
381
382 // Create elastic omega-problem
385 "ELASTIC_PROBLEM_OMEGA_SCHUR");
395
396 // Create elastic tet_stress-problem
399 "ELASTIC_PROBLEM_SPATIAL_DISP_SCHUR");
403 elementVolumeName.c_str());
406 essentialBcElement.c_str());
410
411 {
412 PetscSection section;
413 CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
414 &section);
415 CHKERR DMSetDefaultSection(dmElastic, section);
416 CHKERR DMSetDefaultGlobalSection(dmElastic, section);
417 CHKERR PetscSectionDestroy(&section);
418 }
419
421}
@ QUIET
Definition: definitions.h:208
@ NOISY
Definition: definitions.h:211
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1111
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:221
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:485
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:444
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:118
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:244
MoFEMErrorCode removeDofsOnEntities(const std::string problem_name, const std::string field_name, const Range ents, const int lo_coeff=0, const int hi_coeff=MAX_DOFS_ON_ENTITY, const int lo_order=0, const int hi_order=100, int verb=VERBOSE, const bool debug=false)
Remove DOFs from problem.
auto bit
set bit
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMoFEM.cpp:426
auto createSmartDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
SmartPetscObj< DM > dmElasticSchurBubble
Sub problem of dmElastic Schur.
SmartPetscObj< DM > dmElasticSchurStreach
Sub problem of dmElastic Schur.
SmartPetscObj< DM > dmElasticSchurOmega
Sub problem of dmElastic Schur.
SmartPetscObj< DM > dmElasticSchurSpatialDisp
Sub problem of dmElastic Schur.
SmartPetscObj< DM > dM
Coupled problem all fields.
SmartPetscObj< DM > dmElastic
Elastic problem.
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTraction
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MPI_Comm & get_comm() const =0
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
Problem manager is used to build and partition problems.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

◆ addFields()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::addFields ( const EntityHandle  meshset = 0)
Examples
EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 86 of file EshelbianPlasticity.cpp.

86 {
88
89 Range tets;
90 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
91 Range tets_skin_part;
92 Skinner skin(&mField.get_moab());
93 CHKERR skin.find_skin(0, tets, false, tets_skin_part);
94 ParallelComm *pcomm =
95 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
96 Range tets_skin;
97 CHKERR pcomm->filter_pstatus(tets_skin_part,
98 PSTATUS_SHARED | PSTATUS_MULTISHARED,
99 PSTATUS_NOT, -1, &tets_skin);
100
101 auto subtract_faces_where_displacements_are_applied =
102 [&](const std::string disp_block_set_name) {
105 if (it->getName().compare(0, disp_block_set_name.length(),
106 disp_block_set_name) == 0) {
107 Range faces;
108 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2,
109 faces, true);
110 tets_skin = subtract(tets_skin, faces);
111 }
112 }
114 };
115 CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_DISP_BC");
116 CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_ROTATION_BC");
117 CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_TRACTION_BC");
118
119 Range faces;
120 CHKERR mField.get_moab().get_adjacencies(tets, 2, true, faces,
121 moab::Interface::UNION);
122 Range faces_not_on_the_skin = subtract(faces, tets_skin);
123
124 auto add_hdiv_field = [&](const std::string field_name, const int order,
125 const int dim) {
128 MB_TAG_SPARSE, MF_ZERO);
130 CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
131 CHKERR mField.set_field_order(faces_not_on_the_skin, field_name, order);
134 };
135
136 auto add_hdiv_rt_field = [&](const std::string field_name, const int order,
137 const int dim) {
140 MB_TAG_DENSE, MF_ZERO);
142 CHKERR mField.set_field_order(meshset, MBTET, field_name, 0);
145 };
146
147 auto add_l2_field = [this, meshset](const std::string field_name,
148 const int order, const int dim) {
151 MB_TAG_DENSE, MF_ZERO);
153 CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
155 };
156
157 auto add_h1_field = [this, meshset](const std::string field_name,
158 const int order, const int dim) {
161 MB_TAG_DENSE, MF_ZERO);
163 CHKERR mField.set_field_order(meshset, MBVERTEX, field_name, 1);
164 CHKERR mField.set_field_order(meshset, MBEDGE, field_name, order);
165 CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
166 CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
168 };
169
170 auto add_bubble_field = [this, meshset](const std::string field_name,
171 const int order, const int dim) {
174 MF_ZERO);
175 // Modify field
176 auto field_ptr = mField.get_field_structure(field_name);
177 auto field_order_table =
178 const_cast<Field *>(field_ptr)->getFieldOrderTable();
179 auto get_cgg_bubble_order_zero = [](int p) { return 0; };
180 auto get_cgg_bubble_order_tet = [](int p) {
182 };
183 field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
184 field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
185 field_order_table[MBTRI] = get_cgg_bubble_order_zero;
186 field_order_table[MBTET] = get_cgg_bubble_order_tet;
188 CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
189 CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
191 };
192
193 // spatial fields
194 CHKERR add_hdiv_field(piolaStress, spaceOrder, 3);
195 CHKERR add_bubble_field(bubbleField, spaceOrder, 1);
196 CHKERR add_l2_field(spatialDisp, spaceOrder - 1, 3);
197 CHKERR add_l2_field(rotAxis, spaceOrder - 1, 3);
198 CHKERR add_l2_field(streachTensor, spaceOrder, 6);
199
200 // material fields
201 CHKERR add_hdiv_field(eshelbyStress, materialOrder, 3);
202 CHKERR add_l2_field(materialGradient, materialOrder - 1, 9);
203 // CHKERR add_l2_field(materialDisp, materialOrder - 1, 3);
204 // CHKERR add_l2_field(tauField, materialOrder - 1, 1);
205 // CHKERR add_l2_field(lambdaField, materialOrder - 1, 1);
206
207 // Add history filedes
208 CHKERR add_l2_field(materialGradient + "0", materialOrder - 1, 9);
209
211
213}
static Index< 'p', 3 > p
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
const int dim
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
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.
Provide data structure for (tensor) field approximation.

◆ addMaterial_HMHHStVenantKirchhoff()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::addMaterial_HMHHStVenantKirchhoff ( const int  tape,
const double  lambda,
const double  mu,
const double  sigma_y 
)

Definition at line 633 of file EshelbianADOL-C.cpp.

635 {
637 physicalEquations = boost::shared_ptr<HMHStVenantKirchhoff>(
638 new HMHStVenantKirchhoff(lambda, mu, sigma_y));
639 CHKERR physicalEquations->recordTape(tape, nullptr);
641}
constexpr double lambda
constexpr double mu
boost::shared_ptr< PhysicalEquations > physicalEquations

◆ addMaterial_HMHMooneyRivlin()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::addMaterial_HMHMooneyRivlin ( const int  tape,
const double  alpha,
const double  beta,
const double  lambda,
const double  sigma_y 
)
Examples
ep.cpp.

Definition at line 643 of file EshelbianADOL-C.cpp.

645 {
647 physicalEquations = boost::shared_ptr<HMHPMooneyRivlinWriggersEq63>(
648 new HMHPMooneyRivlinWriggersEq63(alpha, beta, lambda, sigma_y));
649 CHKERR physicalEquations->recordTape(tape, nullptr);
651}

◆ addVolumeFiniteElement()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::addVolumeFiniteElement ( const EntityHandle  meshset = 0)
Examples
EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 216 of file EshelbianPlasticity.cpp.

216 {
218
219 // set finite element fields
220 auto add_field_to_fe = [this](const std::string fe,
221 const std::string field_name) {
227 };
228
233
234 CHKERR add_field_to_fe(elementVolumeName, piolaStress);
235 CHKERR add_field_to_fe(elementVolumeName, bubbleField);
236 CHKERR add_field_to_fe(elementVolumeName, eshelbyStress);
237 CHKERR add_field_to_fe(elementVolumeName, streachTensor);
238 CHKERR add_field_to_fe(elementVolumeName, rotAxis);
239 CHKERR add_field_to_fe(elementVolumeName, spatialDisp);
240 CHKERR add_field_to_fe(elementVolumeName, streachTensor);
241 CHKERR add_field_to_fe(elementVolumeName, materialGradient);
242
244 materialGradient + "0");
245 }
246
247 // build finite elements data structures
249
251}
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.

◆ getBc()

template<typename BC >
MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getBc ( boost::shared_ptr< BC > &  bc_vec_ptr,
const std::string  block_set_name,
const int  nb_attributes 
)
inline

Definition at line 1221 of file EshelbianPlasticity.hpp.

1223 {
1226 auto block_name = it->getName();
1227 if (block_name.compare(0, block_set_name.length(), block_set_name) == 0) {
1228 std::vector<double> block_attributes;
1229 CHKERR it->getAttributes(block_attributes);
1230 if (block_attributes.size() != nb_attributes) {
1231 SETERRQ3(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1232 "In block %s expected %d attributes, but given %d",
1233 it->getName().c_str(), nb_attributes,
1234 block_attributes.size());
1235 }
1236 Range faces;
1237 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
1238 true);
1239 bc_vec_ptr->emplace_back(block_name, block_attributes, faces);
1240 }
1241 }
1243 }
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31

◆ getOptions()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getOptions ( )
Examples
EshelbianPlasticity.cpp.

Definition at line 53 of file EshelbianPlasticity.cpp.

53 {
55 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
56 "none");
57
58 spaceOrder = 1;
59 CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
60 spaceOrder, &spaceOrder, PETSC_NULL);
61 materialOrder = 1;
62 CHKERR PetscOptionsInt("-material_order", "approximation oder for material",
63 "", materialOrder, &materialOrder, PETSC_NULL);
64
65 alphaU = 0;
66 CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alphaU,
67 &alphaU, PETSC_NULL);
68
69 alphaW = 0;
70 CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alphaW,
71 &alphaW, PETSC_NULL);
72
73 alphaRho = 0;
74 CHKERR PetscOptionsScalar("-density_alpha_rho", "density", "", alphaRho,
75 &alphaRho, PETSC_NULL);
76
77 precEps = 1e-3;
78 CHKERR PetscOptionsScalar("-preconditioner_eps", "preconditioner_eps", "",
79 precEps, &precEps, PETSC_NULL);
80
81 ierr = PetscOptionsEnd();
84}
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483

◆ getSpatialDispBc()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getSpatialDispBc ( )
inline
Examples
ep.cpp.

Definition at line 1245 of file EshelbianPlasticity.hpp.

1245 {
1246 bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
1247 return getBc(bcSpatialDispVecPtr, "SPATIAL_DISP_BC", 6);
1248 }
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
MoFEMErrorCode getBc(boost::shared_ptr< BC > &bc_vec_ptr, const std::string block_set_name, const int nb_attributes)

◆ getSpatialRotationBc()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getSpatialRotationBc ( )
inline
Examples
ep.cpp.

Definition at line 1250 of file EshelbianPlasticity.hpp.

1250 {
1251 bcSpatialRotationVecPtr = boost::make_shared<BcRotVec>();
1252 return getBc(bcSpatialRotationVecPtr, "SPATIAL_ROTATION_BC", 4);
1253 }
boost::shared_ptr< BcRotVec > bcSpatialRotationVecPtr

◆ getSpatialTractionBc()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getSpatialTractionBc ( )
inline
Examples
ep.cpp.

Definition at line 1255 of file EshelbianPlasticity.hpp.

1255 {
1256 bcSpatialTraction = boost::make_shared<TractionBcVec>();
1257 return getBc(bcSpatialTraction, "SPATIAL_TRACTION_BC", 6);
1258 }
boost::shared_ptr< TractionBcVec > bcSpatialTraction

◆ getSpatialTractionFreeBc()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getSpatialTractionFreeBc ( const EntityHandle  meshset = 0)
inline
Examples
ep.cpp.

Definition at line 1265 of file EshelbianPlasticity.hpp.

1265 {
1267 boost::shared_ptr<TractionFreeBc>(new TractionFreeBc());
1268 return getTractionFreeBc(meshset, bcSpatialFreeTraction, "SPATIAL_DISP_BC",
1269 "SPATIAL_ROTATION_BC");
1270 }
std::vector< Range > TractionFreeBc
MoFEMErrorCode getTractionFreeBc(const EntityHandle meshset, boost::shared_ptr< TractionFreeBc > &bc_ptr, const std::string disp_block_set_name, const std::string rot_block_set_name)

◆ getTractionFreeBc()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getTractionFreeBc ( const EntityHandle  meshset,
boost::shared_ptr< TractionFreeBc > &  bc_ptr,
const std::string  disp_block_set_name,
const std::string  rot_block_set_name 
)
Examples
EshelbianPlasticity.cpp.

Definition at line 454 of file EshelbianPlasticity.cpp.

457 {
459 Range tets;
460 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
461 Range tets_skin_part;
462 Skinner skin(&mField.get_moab());
463 CHKERR skin.find_skin(0, tets, false, tets_skin_part);
464 ParallelComm *pcomm =
465 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
466 Range tets_skin;
467 CHKERR pcomm->filter_pstatus(tets_skin_part,
468 PSTATUS_SHARED | PSTATUS_MULTISHARED,
469 PSTATUS_NOT, -1, &tets_skin);
470
471 bc_ptr->resize(3);
472 for (int dd = 0; dd != 3; ++dd)
473 (*bc_ptr)[dd] = tets_skin;
474
476 if (it->getName().compare(0, disp_block_set_name.length(),
477 disp_block_set_name) == 0) {
478 std::vector<double> block_attributes;
479 CHKERR it->getAttributes(block_attributes);
480 if (block_attributes.size() != 6) {
481 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
482 "In block %s six attributes are required for given BC "
483 "blockset (3 values + "
484 "3 flags) != %d",
485 it->getName().c_str(), block_attributes.size());
486 }
487 Range faces;
488 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
489 true);
490 if (block_attributes[3] != 0)
491 (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
492 if (block_attributes[4] != 0)
493 (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
494 if (block_attributes[5] != 0)
495 (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
496 }
497 if (it->getName().compare(0, rot_block_set_name.length(),
498 rot_block_set_name) == 0) {
499 Range faces;
500 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
501 true);
502 (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
503 (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
504 (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
505 }
506 }
507
508 // for (int dd = 0; dd != 3; ++dd) {
509 // EntityHandle meshset;
510 // CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
511 // CHKERR mField.get_moab().add_entities(meshset, (*bc_ptr)[dd]);
512 // std::string file_name = disp_block_set_name +
513 // "_traction_free_bc_" + boost::lexical_cast<std::string>(dd) + ".vtk";
514 // CHKERR mField.get_moab().write_file(file_name.c_str(), " VTK ", "",
515 // &meshset, 1);
516 // CHKERR mField.get_moab().delete_entities(&meshset, 1);
517 // }
518
520}
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33

◆ postProcessResults()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::postProcessResults ( const int  tag,
const std::string  file 
)
Examples
EshelbianPlasticity.cpp.

Definition at line 1471 of file EshelbianPlasticity.cpp.

1472 {
1474
1475 if (!dataAtPts) {
1476 dataAtPts =
1477 boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
1478 }
1479
1481 post_proc.getUserPolynomialBase() =
1482 boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1483 post_proc.getOpPtrVector().push_back(new OpL2Transform());
1484
1485 CHKERR post_proc.generateReferenceElementMesh();
1486 post_proc.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
1487 piolaStress, dataAtPts->getApproxPAtPts()));
1488 post_proc.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
1489 bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
1490 post_proc.getOpPtrVector().push_back(
1492 streachTensor, dataAtPts->getLogStreachTensorAtPts(), MBTET));
1493 post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1494 rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
1495 post_proc.getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
1496 materialGradient, dataAtPts->getBigGAtPts(), MBTET));
1497 post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1498 spatialDisp, dataAtPts->getSmallWAtPts(), MBTET));
1499
1500 // evaluate derived quantities
1501 post_proc.getOpPtrVector().push_back(
1502 new OpCalculateRotationAndSpatialGradient(rotAxis, dataAtPts));
1503
1504 // evaluate integration points
1505 post_proc.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
1506 spatialDisp, tag, true, false, dataAtPts, physicalEquations));
1507 post_proc.getOpPtrVector().push_back(new OpPostProcDataStructure(
1508 post_proc.postProcMesh, post_proc.mapGaussPts, spatialDisp, dataAtPts));
1509
1511 CHKERR post_proc.writeFile(file.c_str());
1513}
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:574
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
Calculate symmetric tensor field values at integration pts.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Post processing.

◆ query_interface()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::query_interface ( boost::typeindex::type_index  type_index,
UnknownInterface **  iface 
) const

Getting interface of core database.

Parameters
uuidunique ID of interface
ifacereturned pointer to interface
Returns
error code
Examples
EshelbianPlasticity.cpp.

Definition at line 20 of file EshelbianPlasticity.cpp.

21 {
22 *iface = const_cast<EshelbianCore *>(this);
23 return 0;
24}
EshelbianCore(MoFEM::Interface &m_field)

◆ setBaseVolumeElementOps()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::setBaseVolumeElementOps ( const int  tag,
const bool  do_rhs,
const bool  do_lhs,
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > &  fe 
)
Examples
EshelbianPlasticity.cpp.

Definition at line 630 of file EshelbianPlasticity.cpp.

632 {
634 fe = boost::make_shared<EpElement<VolumeElementForcesAndSourcesCore>>(mField);
635
636 fe->getUserPolynomialBase() =
637 boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
638 fe->getOpPtrVector().push_back(new OpL2Transform());
639
640 // set integration rule
641 fe->getRuleHook = VolRule();
642
643 if (!dataAtPts) {
644 dataAtPts =
645 boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
646 dataAtPts->physicsPtr = physicalEquations;
647 }
648
649 // calculate fields values
650 fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
651 piolaStress, dataAtPts->getApproxPAtPts()));
652 fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
653 bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
654 fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
655 piolaStress, dataAtPts->getDivPAtPts()));
656 fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
657 eshelbyStress, dataAtPts->getApproxSigmaAtPts()));
658 fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
659 eshelbyStress, dataAtPts->getDivSigmaAtPts()));
660 fe->getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
661 streachTensor, dataAtPts->getLogStreachTensorAtPts(), MBTET));
662 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
663 rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
664 fe->getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
665 materialGradient, dataAtPts->getBigGAtPts(), MBTET));
666 fe->getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
667 materialGradient + "0", dataAtPts->getBigG0AtPts(), MBTET));
668 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
669 spatialDisp, dataAtPts->getSmallWAtPts(), MBTET));
670
671 // velocities
672 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
673 spatialDisp, dataAtPts->getSmallWDotAtPts(), MBTET));
674 fe->getOpPtrVector().push_back(
676 streachTensor, dataAtPts->getLogStreachDotTensorAtPts(), MBTET));
677 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
678 rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
679
680 // acceleration
681 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
682 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDotDot<3>(
683 spatialDisp, dataAtPts->getSmallWDotDotAtPts(), MBTET));
684 }
685
686 // calculate other derived quantities
687 fe->getOpPtrVector().push_back(
688 new OpCalculateRotationAndSpatialGradient(rotAxis, dataAtPts));
689
690 // evaluate integration points
691 fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
692 spatialDisp, tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
693
695}
Calculate divergence of tonsorial field using vectorial base.
Calculate symmetric tensor field rates ant integratio pts.
Approximate field values for given petsc vector.
Set integration rule.

◆ setElasticElementOps()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::setElasticElementOps ( const int  tag)
Examples
EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 811 of file EshelbianPlasticity.cpp.

811 {
817}
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > elasticFeLhs
MoFEMErrorCode setGenericFaceElementOps(const bool add_elastic, const bool add_material, boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > &fe_rhs, boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > &fe_lhs)
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > elasticBcRhs
MoFEMErrorCode setGenericVolumeElementOps(const int tag, const bool add_elastic, const bool add_material, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > &fe_rhs, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > &fe_lhs)
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > elasticFeRhs
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > elasticBcLhs

◆ setElasticElementToTs()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::setElasticElementToTs ( DM  dm)
Examples
EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 819 of file EshelbianPlasticity.cpp.

819 {
821 boost::shared_ptr<FEMethod> null;
822 boost::shared_ptr<EpElement<FeTractionBc>> spatial_traction_bc(
823 new EpElement<FeTractionBc>(mField, piolaStress, bcSpatialTraction));
824 schurAssembly = boost::make_shared<EpFEMethod>();
825
826 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
827 CHKERR DMMoFEMTSSetI2Function(dm, DM_NO_ELEMENT, null, spatial_traction_bc,
828 null);
830 null);
832 null);
834 spatial_traction_bc);
837 null);
839 null);
841 } else {
842 CHKERR DMMoFEMTSSetIFunction(dm, DM_NO_ELEMENT, null, spatial_traction_bc,
843 null);
845 null);
847 null);
849 spatial_traction_bc);
852 null);
854 null);
856 }
858}
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMoFEM.cpp:788
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMoFEM.cpp:841
PetscErrorCode DMMoFEMTSSetI2Jacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMoFEM.cpp:1005
PetscErrorCode DMMoFEMTSSetI2Function(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS implicit function evaluation function
Definition: DMMoFEM.cpp:963
boost::shared_ptr< EpFEMethod > schurAssembly

◆ setGenericFaceElementOps()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::setGenericFaceElementOps ( const bool  add_elastic,
const bool  add_material,
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > &  fe_rhs,
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > &  fe_lhs 
)
Examples
EshelbianPlasticity.cpp.

Definition at line 781 of file EshelbianPlasticity.cpp.

784 {
786
787 fe_rhs =
788 boost::make_shared<EpElement<FaceElementForcesAndSourcesCore>>(mField);
789 fe_lhs =
790 boost::make_shared<EpElement<FaceElementForcesAndSourcesCore>>(mField);
791
792 // set integration rule
793 fe_rhs->getRuleHook = FaceRule();
794 fe_lhs->getRuleHook = FaceRule();
795
796 fe_rhs->getOpPtrVector().push_back(
798 fe_lhs->getOpPtrVector().push_back(
800
801 if (add_elastic) {
802 fe_rhs->getOpPtrVector().push_back(
804 fe_rhs->getOpPtrVector().push_back(
805 new OpRotationBc(piolaStress, dataAtPts, bcSpatialRotationVecPtr));
806 }
807
809}
Set integration rule to boundary elements.
transform Hdiv base fluxes from reference element to physical triangle

◆ setGenericVolumeElementOps()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::setGenericVolumeElementOps ( const int  tag,
const bool  add_elastic,
const bool  add_material,
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > &  fe_rhs,
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > &  fe_lhs 
)
Examples
EshelbianPlasticity.cpp.

Definition at line 697 of file EshelbianPlasticity.cpp.

700 {
702
703 // Right hand side
704 CHKERR setBaseVolumeElementOps(tag, true, false, fe_rhs);
705
706 // elastic
707 if (add_elastic) {
708 fe_rhs->getOpPtrVector().push_back(
709 new OpSpatialEquilibrium(spatialDisp, dataAtPts, alphaW, alphaRho));
710 fe_rhs->getOpPtrVector().push_back(
711 new OpSpatialRotation(rotAxis, dataAtPts));
712 fe_rhs->getOpPtrVector().push_back(
713 new OpSpatialPhysical(streachTensor, dataAtPts, alphaU));
714 fe_rhs->getOpPtrVector().push_back(
715 new OpSpatialConsistencyP(piolaStress, dataAtPts));
716 fe_rhs->getOpPtrVector().push_back(
717 new OpSpatialConsistencyBubble(bubbleField, dataAtPts));
718 fe_rhs->getOpPtrVector().push_back(
719 new OpSpatialConsistencyDivTerm(piolaStress, dataAtPts));
720 }
721
722 // Left hand side
723 CHKERR setBaseVolumeElementOps(tag, true, true, fe_lhs);
724
725 // elastic
726 if (add_elastic) {
727
728 // Schur
729 fe_lhs->getOpPtrVector().push_back(
730 new OpSpatialSchurBegin(spatialDisp, dataAtPts));
731
732 fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_du(
734 fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
736
737 fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
739 fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
741
742 fe_lhs->getOpPtrVector().push_back(
743 new OpSpatialConsistency_dP_domega(piolaStress, rotAxis, dataAtPts));
744 fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
746
747 fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
749 fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
751
752 fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
753 rotAxis, piolaStress, dataAtPts, false));
754 fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
755 rotAxis, bubbleField, dataAtPts, false));
756 fe_lhs->getOpPtrVector().push_back(
757 new OpSpatialRotation_domega_domega(rotAxis, rotAxis, dataAtPts));
758
759 // Schur
760 dataAtPts->ooMatPtr = boost::make_shared<MatrixDouble>();
761 fe_lhs->getOpPtrVector().push_back(
762 new OpSpatialPreconditionMass(rotAxis, dataAtPts->ooMatPtr));
763 if (alphaW < std::numeric_limits<double>::epsilon() &&
764 alphaRho < std::numeric_limits<double>::epsilon()) {
765 dataAtPts->wwMatPtr = boost::make_shared<MatrixDouble>();
766 fe_lhs->getOpPtrVector().push_back(
767 new OpSpatialPreconditionMass(spatialDisp, dataAtPts->wwMatPtr));
768 } else {
769 dataAtPts->wwMatPtr.reset();
770 }
771 fe_lhs->getOpPtrVector().push_back(
772 new OpSpatialSchurEnd(spatialDisp, dataAtPts, precEps));
773
774 if (add_material) {
775 }
776 }
777
779}
MoFEMErrorCode setBaseVolumeElementOps(const int tag, const bool do_rhs, const bool do_lhs, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > &fe)

◆ setUpTSElastic()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::setUpTSElastic ( TS  ts,
Mat  m,
Vec  f 
)
Examples
EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 860 of file EshelbianPlasticity.cpp.

860 {
862 boost::shared_ptr<TsCtx> ts_ctx;
864
865 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
866 CHKERR TSSetI2Function(ts, f, PETSC_NULL, PETSC_NULL);
867 CHKERR TSSetI2Jacobian(ts, m, m, PETSC_NULL, PETSC_NULL);
868 } else {
869 CHKERR TSSetIFunction(ts, f, PETSC_NULL, PETSC_NULL);
870 CHKERR TSSetIJacobian(ts, m, m, PETSC_NULL, PETSC_NULL);
871 }
872
873 CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx.get(), PETSC_NULL);
874
875 auto add_schur_streach_op = [this](auto &list, SmartPetscObj<Mat> S,
876 SmartPetscObj<AO> aoS) {
878 for (auto &fe : list)
879 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
880 fe_cast->addStreachSchurMatrix(S, aoS);
881 else
882 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
884 };
885
886 auto add_schur_streach_pre = [this](auto &list, SmartPetscObj<Mat> S,
887 SmartPetscObj<AO> aoS) {
889 for (auto &fe : list)
890 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
891 fe_cast->addStreachSchurMatrix(S, aoS);
892 else
893 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
895 };
896
897 auto add_schur_bubble_op = [this](auto &list, SmartPetscObj<Mat> S,
898 SmartPetscObj<AO> aoS) {
900 for (auto &fe : list)
901 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
902 fe_cast->addBubbleSchurMatrix(S, aoS);
903 else
904 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
906 };
907
908 auto add_schur_bubble_pre = [this](auto &list, SmartPetscObj<Mat> S,
909 SmartPetscObj<AO> aoS) {
911 for (auto &fe : list)
912 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
913 fe_cast->addBubbleSchurMatrix(S, aoS);
914 else
915 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
917 };
918
919 auto add_schur_spatial_disp_op = [this](auto &list, SmartPetscObj<Mat> S,
920 SmartPetscObj<AO> aoS) {
922 for (auto &fe : list)
923 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
924 fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
925 else
926 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
928 };
929
930 auto add_schur_spatial_disp_pre = [this](auto &list, SmartPetscObj<Mat> S,
931 SmartPetscObj<AO> aoS) {
933 for (auto &fe : list)
934 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
935 fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
936 else
937 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
939 };
940
941 auto add_schur_omega_op = [this](auto &list, SmartPetscObj<Mat> S,
942 SmartPetscObj<AO> aoS) {
944 for (auto &fe : list)
945 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
946 fe_cast->addOmegaSchurMatrix(S, aoS);
947 else
948 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
950 };
951
952 auto add_schur_omega_pre = [this](auto &list, SmartPetscObj<Mat> S,
955 for (auto &fe : list)
956 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
957 fe_cast->addOmegaSchurMatrix(S, ao);
958 else
959 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
961 };
962
963 const MoFEM::Problem *schur_streach_prb_ptr;
964 CHKERR DMMoFEMGetProblemPtr(dmElasticSchurStreach, &schur_streach_prb_ptr);
965 if (auto sub_data = schur_streach_prb_ptr->subProblemData) {
968 SmartPetscObj<AO> aoSuu = sub_data->getSmartRowMap();
969
970 CHKERR add_schur_streach_op(ts_ctx->loopsIJacobian, Suu, aoSuu);
971 CHKERR add_schur_streach_pre(ts_ctx->preProcessIJacobian, Suu, aoSuu);
972 CHKERR add_schur_streach_pre(ts_ctx->postProcessIJacobian, Suu, aoSuu);
973
974 const MoFEM::Problem *schur_bubble_prb_ptr;
975 CHKERR DMMoFEMGetProblemPtr(dmElasticSchurBubble, &schur_bubble_prb_ptr);
976 if (auto bubble_data = schur_bubble_prb_ptr->subProblemData) {
977 SmartPetscObj<Mat> SBubble;
979 SmartPetscObj<AO> aoSBubble = bubble_data->getSmartRowMap();
980 CHKERR add_schur_bubble_op(ts_ctx->loopsIJacobian, SBubble,
981 aoSBubble);
982 CHKERR add_schur_bubble_pre(ts_ctx->preProcessIJacobian, SBubble,
983 aoSBubble);
984 CHKERR add_schur_bubble_pre(ts_ctx->postProcessIJacobian, SBubble,
985 aoSBubble);
986
987 const MoFEM::Problem *schur_omega_prb_ptr;
988 CHKERR DMMoFEMGetProblemPtr(dmElasticSchurOmega, &schur_omega_prb_ptr);
989 if (auto tet_stress_data = schur_omega_prb_ptr->subProblemData) {
990 SmartPetscObj<Mat> SOmega;
992 SmartPetscObj<AO> aoSOmega = tet_stress_data->getSmartRowMap();
993 CHKERR add_schur_omega_op(ts_ctx->loopsIJacobian, SOmega,
994 aoSOmega);
995 CHKERR add_schur_omega_pre(ts_ctx->preProcessIJacobian, SOmega,
996 aoSOmega);
997 CHKERR add_schur_omega_pre(ts_ctx->postProcessIJacobian, SOmega,
998 aoSOmega);
999
1000 const MoFEM::Problem *schur_spatial_disp_prb_ptr;
1002 &schur_spatial_disp_prb_ptr);
1003 if (auto spatial_disp_data =
1004 schur_spatial_disp_prb_ptr->subProblemData) {
1007 SmartPetscObj<AO> aoSw = spatial_disp_data->getSmartRowMap();
1008 CHKERR add_schur_spatial_disp_op(ts_ctx->loopsIJacobian, Sw,
1009 aoSw);
1010 CHKERR add_schur_spatial_disp_pre(ts_ctx->preProcessIJacobian, Sw,
1011 aoSw);
1012 CHKERR add_schur_spatial_disp_pre(ts_ctx->postProcessIJacobian, Sw,
1013 aoSw);
1014 } else
1016 "Problem does not have sub-problem data");
1017
1018 } else
1020 "Problem does not have sub-problem data");
1021
1022 } else
1024 "Problem does not have sub-problem data");
1025
1026 } else
1028 "Problem does not have sub-problem data");
1029
1030 struct Monitor : public FEMethod {
1031
1032 using Ele = ForcesAndSourcesCore;
1034 using VolOp = VolumeElementForcesAndSourcesCore::UserDataOperator;
1035 using SetPtsData = FieldEvaluatorInterface::SetPtsData;
1036
1037 EshelbianCore &eP;
1038 boost::shared_ptr<SetPtsData> dataFieldEval;
1039 boost::shared_ptr<VolEle> volPostProcEnergy;
1040 boost::shared_ptr<double> gEnergy;
1041
1042 Monitor(EshelbianCore &ep)
1043 : eP(ep),
1044 dataFieldEval(ep.mField.getInterface<FieldEvaluatorInterface>()
1045 ->getData<VolEle>()),
1046 volPostProcEnergy(new VolEle(ep.mField)), gEnergy(new double) {
1047 ierr = ep.mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
1048 dataFieldEval, "EP");
1049 CHKERRABORT(PETSC_COMM_WORLD, ierr);
1050
1051 auto no_rule = [](int, int, int) { return -1; };
1052
1053 auto set_element_for_field_eval = [&]() {
1054 boost::shared_ptr<Ele> vol_ele(dataFieldEval->feMethodPtr.lock());
1055 vol_ele->getRuleHook = no_rule;
1056 vol_ele->getUserPolynomialBase() =
1057 boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1058 vol_ele->getOpPtrVector().push_back(new OpL2Transform());
1059
1060 vol_ele->getOpPtrVector().push_back(
1062 eP.piolaStress, eP.dataAtPts->getApproxPAtPts()));
1063 vol_ele->getOpPtrVector().push_back(
1065 eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
1066 vol_ele->getOpPtrVector().push_back(
1068 eP.streachTensor, eP.dataAtPts->getLogStreachTensorAtPts(),
1069 MBTET));
1070 vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1071 eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
1072 vol_ele->getOpPtrVector().push_back(
1074 eP.materialGradient, eP.dataAtPts->getBigGAtPts(), MBTET));
1075 vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1076 eP.spatialDisp, eP.dataAtPts->getSmallWAtPts(), MBTET));
1077 vol_ele->getOpPtrVector().push_back(
1078 new OpCalculateRotationAndSpatialGradient(eP.rotAxis,
1079 eP.dataAtPts));
1080 };
1081
1082 auto set_element_for_post_process = [&]() {
1083 volPostProcEnergy->getRuleHook = VolRule();
1084 volPostProcEnergy->getUserPolynomialBase() =
1085 boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1086 volPostProcEnergy->getOpPtrVector().push_back(new OpL2Transform());
1087
1088 volPostProcEnergy->getOpPtrVector().push_back(
1090 eP.piolaStress, eP.dataAtPts->getApproxPAtPts()));
1091 volPostProcEnergy->getOpPtrVector().push_back(
1093 eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
1094 volPostProcEnergy->getOpPtrVector().push_back(
1096 eP.streachTensor, eP.dataAtPts->getLogStreachTensorAtPts(),
1097 MBTET));
1098 volPostProcEnergy->getOpPtrVector().push_back(
1100 eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
1101 volPostProcEnergy->getOpPtrVector().push_back(
1103 eP.materialGradient, eP.dataAtPts->getBigGAtPts(), MBTET));
1104 volPostProcEnergy->getOpPtrVector().push_back(
1106 eP.spatialDisp, eP.dataAtPts->getSmallWAtPts(), MBTET));
1107 volPostProcEnergy->getOpPtrVector().push_back(
1108 new OpCalculateRotationAndSpatialGradient(eP.rotAxis,
1109 eP.dataAtPts));
1110 volPostProcEnergy->getOpPtrVector().push_back(
1111 new OpCalculateStrainEnergy(eP.spatialDisp, eP.dataAtPts, gEnergy));
1112 };
1113
1114 set_element_for_field_eval();
1115 set_element_for_post_process();
1116 }
1117
1118 MoFEMErrorCode preProcess() { return 0; }
1119
1120 MoFEMErrorCode operator()() { return 0; }
1121
1122 MoFEMErrorCode postProcess() {
1124
1125 // auto get_str_time = [](auto ts_t) {
1126 // std::ostringstream ss;
1127 // ss << boost::str(boost::format("%d") %
1128 // static_cast<int>(std::ceil(ts_t * 1e6)));
1129 // std::string s = ss.str();
1130 // return s;
1131 // };
1132
1133 auto get_step = [](auto ts_step) {
1134 std::ostringstream ss;
1135 ss << boost::str(boost::format("%d") % static_cast<int>(ts_step));
1136 std::string s = ss.str();
1137 return s;
1138 };
1139
1140 PetscViewer viewer;
1141 CHKERR PetscViewerBinaryOpen(
1142 PETSC_COMM_WORLD, ("restart_" + get_step(ts_step) + ".dat").c_str(),
1143 FILE_MODE_WRITE, &viewer);
1144 CHKERR VecView(ts_u, viewer);
1145 CHKERR PetscViewerDestroy(&viewer);
1146
1147 CHKERR eP.postProcessResults(1, "out_sol_elastic_" + get_step(ts_step) +
1148 ".h5m");
1149
1150 // Loop boundary elements with traction boundary conditions
1151 *gEnergy = 0;
1152 CHKERR eP.mField.loop_finite_elements(problemPtr->getName(), "EP",
1153 *volPostProcEnergy);
1154
1155 double body_energy;
1156 MPI_Allreduce(gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
1157 eP.mField.get_comm());
1158 MOFEM_LOG_C("EP", Sev::inform, "Step %d time %3.4g strain energy %3.6e",
1159 ts_step, ts_t, body_energy);
1160
1161 auto post_proc_at_points = [&](std::array<double, 3> point,
1162 std::string str) {
1164
1165 dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
1166
1167 struct OpPrint : public VolOp {
1168
1169 EshelbianCore &eP;
1170 std::array<double, 3> point;
1171 std::string str;
1172
1173 OpPrint(EshelbianCore &ep, std::array<double, 3> &point,
1174 std::string &str)
1175 : VolOp(ep.spatialDisp, VolOp::OPROW), eP(ep), point(point),
1176 str(str) {}
1177
1178 MoFEMErrorCode doWork(int side, EntityType type,
1179 DataForcesAndSourcesCore::EntData &data) {
1181 if (type == MBTET) {
1182 if (getGaussPts().size2()) {
1183
1184 auto t_h = getFTensor2FromMat<3, 3>(eP.dataAtPts->hAtPts);
1185 auto t_approx_P =
1186 getFTensor2FromMat<3, 3>(eP.dataAtPts->approxPAtPts);
1187
1191 const double jac = determinantTensor3by3(t_h);
1193 t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
1194
1195 auto add = [&]() {
1196 std::ostringstream s;
1197 s << str << " elem " << getFEEntityHandle() << " ";
1198 return s.str();
1199 };
1200
1201 auto print_tensor = [](auto &t) {
1202 std::ostringstream s;
1203 s << t;
1204 return s.str();
1205 };
1206
1207 std::ostringstream print;
1208 MOFEM_LOG("EPSYNC", Sev::inform)
1209 << add() << "comm rank " << eP.mField.get_comm_rank();
1210 MOFEM_LOG("EPSYNC", Sev::inform)
1211 << add() << "point " << getVectorAdaptor(point.data(), 3);
1212 MOFEM_LOG("EPSYNC", Sev::inform)
1213 << add() << "coords at gauss pts " << getCoordsAtGaussPts();
1214 MOFEM_LOG("EPSYNC", Sev::inform)
1215 << add() << "w " << eP.dataAtPts->wAtPts;
1216 MOFEM_LOG("EPSYNC", Sev::inform)
1217 << add() << "Piola " << eP.dataAtPts->approxPAtPts;
1218 MOFEM_LOG("EPSYNC", Sev::inform)
1219 << add() << "Cauchy " << print_tensor(t_cauchy);
1220 }
1221 }
1223 }
1224 };
1225
1226 if (auto fe_ptr = dataFieldEval->feMethodPtr.lock()) {
1227
1228 fe_ptr->getOpPtrVector().push_back(new OpPrint(eP, point, str));
1229 CHKERR eP.mField.getInterface<FieldEvaluatorInterface>()
1230 ->evalFEAtThePoint3D(
1231 point.data(), 1e-12, problemPtr->getName(), "EP",
1232 dataFieldEval, eP.mField.get_comm_rank(),
1233 eP.mField.get_comm_rank(), nullptr, MF_EXIST, QUIET);
1234 fe_ptr->getOpPtrVector().pop_back();
1235 }
1236
1238 };
1239
1240 // Points for Cook beam
1241 std::array<double, 3> pointA = {48., 60., 5.};
1242 CHKERR post_proc_at_points(pointA, "Point A");
1243 MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
1244
1245 std::array<double, 3> pointB = {48. / 2., 44. + (60. - 44.) / 2., 0.};
1246 CHKERR post_proc_at_points(pointB, "Point B");
1247 MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
1248
1249 std::array<double, 3> pointC = {48. / 2., (44. - 0.) / 2., 0.};
1250 CHKERR post_proc_at_points(pointC, "Point C");
1251 MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
1252
1254 }
1255 };
1256
1257 boost::shared_ptr<FEMethod> monitor_ptr(new Monitor(*this));
1258 ts_ctx->getLoopsMonitor().push_back(
1260
1261 CHKERR TSAppendOptionsPrefix(ts, "elastic_");
1262 CHKERR TSSetFromOptions(ts);
1264}
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
@ MF_EXIST
Definition: definitions.h:100
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:414
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1185
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMoFEM.cpp:1130
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
FTensor::Index< 'i', SPACE_DIM > i
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1876
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
double f(const double v)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:31
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:219
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1352
constexpr double t
plate stiffness
Definition: plate.cpp:59
structure for User Loop Methods on finite elements
Field evaluator interface.
structure to get information form mofem into EntitiesFieldData
keeps basic data about problem
boost::shared_ptr< SubProblemData > subProblemData
intrusive_ptr for managing petsc objects
FEMethodsSequence loopsIJacobian
Definition: TsCtx.hpp:27
FEMethodsSequence & getLoopsMonitor()
Get the loops to do Monitor object.
Definition: TsCtx.hpp:108
BasicMethodsSequence preProcessIJacobian
Definition: TsCtx.hpp:32
BasicMethodsSequence postProcessIJacobian
Definition: TsCtx.hpp:33
Monitor solution.

◆ solveElastic()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::solveElastic ( TS  ts,
Vec  x 
)
Examples
EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 1266 of file EshelbianPlasticity.cpp.

1266 {
1268
1269 CHKERR TSSetDM(ts, dmElastic);
1270
1271 SNES snes;
1272 CHKERR TSGetSNES(ts, &snes);
1273
1274 PetscViewerAndFormat *vf;
1275 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
1276 PETSC_VIEWER_DEFAULT, &vf);
1277 CHKERR SNESMonitorSet(
1278 snes,
1279 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
1280 vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
1281
1282 PetscSection section;
1283 CHKERR DMGetDefaultSection(dmElastic, &section);
1284 int num_fields;
1285 CHKERR PetscSectionGetNumFields(section, &num_fields);
1286 for (int ff = 0; ff != num_fields; ff++) {
1287 const char *field_name;
1288 CHKERR PetscSectionGetFieldName(section, ff, &field_name);
1289 MOFEM_LOG_C("EP", Sev::inform, "Field %d name %s", ff, field_name);
1290 }
1291
1292 CHKERR DMoFEMMeshToLocalVector(dmElastic, x, INSERT_VALUES, SCATTER_FORWARD);
1293
1294 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
1295 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
1296
1297
1298
1299 // Adding field split solver
1300
1301 KSP ksp;
1302 CHKERR SNESGetKSP(snes, &ksp);
1303 PC pc;
1304 CHKERR KSPGetPC(ksp, &pc);
1305 PetscBool is_uu_field_split;
1306 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_uu_field_split);
1307 if (is_uu_field_split) {
1308
1309 const MoFEM::Problem *schur_uu_ptr;
1311 if (auto uu_data = schur_uu_ptr->subProblemData) {
1312
1313 const MoFEM::Problem *prb_ptr;
1315 map<std::string, IS> is_map;
1316 for (int ff = 0; ff != num_fields; ff++) {
1317 const char *field_name;
1318 CHKERR PetscSectionGetFieldName(section, ff, &field_name);
1319 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1320 prb_ptr->getName(), ROW, field_name, 0, MAX_DOFS_ON_ENTITY,
1321 &is_map[field_name]);
1322 }
1323 // CHKERR mField.getInterface<ISManager>()
1324 // ->isCreateProblemFieldAndEntityType(
1325 // prb_ptr->getName(), ROW, piolaStress, MBTET, MBTET, 0,
1326 // MAX_DOFS_ON_ENTITY, &is_map["T_STRESS_ON_TETS"]);
1327 // CHKERR mField.getInterface<ISManager>()
1328 // ->isCreateProblemFieldAndEntityType(
1329 // prb_ptr->getName(), ROW, piolaStress, MBTRI, MBTRI, 0,
1330 // MAX_DOFS_ON_ENTITY, &is_map["T_STRESS_ON_TRIS"]);
1331
1332 CHKERR uu_data->getRowIs(&is_map["E_IS_SUU"]);
1333
1334 CHKERR PCFieldSplitSetIS(pc, NULL, is_map[streachTensor]);
1335 CHKERR PCFieldSplitSetIS(pc, NULL, is_map["E_IS_SUU"]);
1336
1337 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER,
1338 schurAssembly->Suu);
1339
1340 CHKERR PCSetUp(pc);
1341 PetscInt n;
1342 KSP *uu_ksp;
1343 CHKERR PCFieldSplitGetSubKSP(pc, &n, &uu_ksp);
1344 PC bubble_pc;
1345 CHKERR KSPGetPC(uu_ksp[1], &bubble_pc);
1346 PetscBool is_bubble_field_split;
1347 PetscObjectTypeCompare((PetscObject)bubble_pc, PCFIELDSPLIT,
1348 &is_bubble_field_split);
1349 if (is_bubble_field_split) {
1350
1351 const MoFEM::Problem *schur_bubble_ptr;
1353 if (auto bubble_data = schur_bubble_ptr->subProblemData) {
1354
1355 CHKERR bubble_data->getRowIs(&is_map["E_IS_BUBBLE"]);
1356
1357 AO uu_ao;
1358 CHKERR uu_data->getRowMap(&uu_ao);
1359
1360 CHKERR AOApplicationToPetscIS(uu_ao, is_map[bubbleField]);
1361 CHKERR PCFieldSplitSetIS(bubble_pc, NULL, is_map[bubbleField]);
1362 CHKERR PCFieldSplitSetIS(bubble_pc, NULL, is_map["E_IS_BUBBLE"]);
1363 CHKERR PCFieldSplitSetSchurPre(
1364 bubble_pc, PC_FIELDSPLIT_SCHUR_PRE_USER, schurAssembly->SBubble);
1365
1366 CHKERR PCSetUp(bubble_pc);
1367 PetscInt bubble_n;
1368 KSP *bubble_ksp;
1369 CHKERR PCFieldSplitGetSubKSP(bubble_pc, &bubble_n, &bubble_ksp);
1370 PC omega_pc;
1371 CHKERR KSPGetPC(bubble_ksp[1], &omega_pc);
1372 PetscBool is_omega_field_split;
1373 PetscObjectTypeCompare((PetscObject)omega_pc, PCFIELDSPLIT,
1374 &is_omega_field_split);
1375
1376 if (is_omega_field_split) {
1377
1378 const MoFEM::Problem *schur_omega_ptr;
1380 if (auto omega_data = schur_omega_ptr->subProblemData) {
1381
1382 AO bubble_ao;
1383 CHKERR bubble_data->getRowMap(&bubble_ao);
1384
1385 CHKERR AOApplicationToPetscIS(uu_ao, is_map[rotAxis]);
1386 CHKERR AOApplicationToPetscIS(bubble_ao, is_map[rotAxis]);
1387 CHKERR omega_data->getRowIs(&is_map["E_IS_OMEGA"]);
1388
1389 CHKERR PCFieldSplitSetIS(omega_pc, NULL, is_map[rotAxis]);
1390 CHKERR PCFieldSplitSetIS(omega_pc, NULL, is_map["E_IS_OMEGA"]);
1391 CHKERR PCFieldSplitSetSchurPre(omega_pc,
1392 PC_FIELDSPLIT_SCHUR_PRE_USER,
1393 schurAssembly->SOmega);
1394
1395 CHKERR PCSetUp(omega_pc);
1396 PetscInt omega_n;
1397 KSP *omega_ksp;
1398 CHKERR PCFieldSplitGetSubKSP(omega_pc, &omega_n, &omega_ksp);
1399 PC w_pc;
1400 CHKERR KSPGetPC(omega_ksp[1], &w_pc);
1401 PetscBool is_w_field_split;
1402 PetscObjectTypeCompare((PetscObject)w_pc, PCFIELDSPLIT,
1403 &is_w_field_split);
1404 if (is_w_field_split) {
1405
1406 const MoFEM::Problem *schur_w_ptr;
1408 &schur_w_ptr);
1409 if (auto w_data = schur_w_ptr->subProblemData) {
1410
1411 AO omega_ao;
1412 CHKERR omega_data->getRowMap(&omega_ao);
1413
1414 CHKERR AOApplicationToPetscIS(uu_ao, is_map[spatialDisp]);
1415 CHKERR AOApplicationToPetscIS(bubble_ao, is_map[spatialDisp]);
1416 CHKERR AOApplicationToPetscIS(omega_ao, is_map[spatialDisp]);
1417 CHKERR w_data->getRowIs(&is_map["E_IS_W"]);
1418
1419 CHKERR PCFieldSplitSetIS(w_pc, NULL, is_map[spatialDisp]);
1420 CHKERR PCFieldSplitSetIS(w_pc, NULL, is_map["E_IS_W"]);
1421 CHKERR PCFieldSplitSetSchurPre(
1422 w_pc, PC_FIELDSPLIT_SCHUR_PRE_USER, schurAssembly->Sw);
1423
1424 CHKERR AODestroy(&omega_ao);
1425 }
1426 }
1427
1428 CHKERR PetscFree(omega_ksp);
1429 }
1430 }
1431 CHKERR PetscFree(bubble_ksp);
1432 CHKERR AODestroy(&uu_ao);
1433 }
1434 }
1435 CHKERR PetscFree(uu_ksp);
1436
1437 for (auto &m : is_map)
1438 CHKERR ISDestroy(&m.second);
1439 }
1440 }
1441
1442 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
1443 Vec xx;
1444 CHKERR VecDuplicate(x, &xx);
1445 CHKERR VecZeroEntries(xx);
1446 CHKERR TS2SetSolution(ts, x, xx);
1447 CHKERR VecDestroy(&xx);
1448 } else {
1449 CHKERR TSSetSolution(ts, x);
1450 }
1451
1452 CHKERR TSSolve(ts, PETSC_NULL);
1453
1454 // CHKERR TSGetSNES(ts, &snes);
1455 int lin_solver_iterations;
1456 CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
1457 MOFEM_LOG("EP", Sev::inform)
1458 << "Number of linear solver iterations " << lin_solver_iterations;
1459
1460 PetscBool test_cook_flg = PETSC_FALSE;
1461 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_cook", &test_cook_flg,
1462 PETSC_NULL);
1463 if (test_cook_flg)
1464 if (lin_solver_iterations != 2)
1465 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1466 "Expected number of iterations is different than expected");
1467
1469}
@ ROW
Definition: definitions.h:123
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
FTensor::Index< 'n', SPACE_DIM > n
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:511
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)

Member Data Documentation

◆ alphaRho

double EshelbianPlasticity::EshelbianCore::alphaRho
Examples
EshelbianPlasticity.cpp.

Definition at line 1210 of file EshelbianPlasticity.hpp.

◆ alphaU

double EshelbianPlasticity::EshelbianCore::alphaU
Examples
EshelbianPlasticity.cpp.

Definition at line 1208 of file EshelbianPlasticity.hpp.

◆ alphaW

double EshelbianPlasticity::EshelbianCore::alphaW
Examples
EshelbianPlasticity.cpp.

Definition at line 1209 of file EshelbianPlasticity.hpp.

◆ bcSpatialDispVecPtr

boost::shared_ptr<BcDispVec> EshelbianPlasticity::EshelbianCore::bcSpatialDispVecPtr
Examples
EshelbianPlasticity.cpp.

Definition at line 1215 of file EshelbianPlasticity.hpp.

◆ bcSpatialFreeTraction

boost::shared_ptr<TractionFreeBc> EshelbianPlasticity::EshelbianCore::bcSpatialFreeTraction
Examples
EshelbianPlasticity.cpp.

Definition at line 1218 of file EshelbianPlasticity.hpp.

◆ bcSpatialRotationVecPtr

boost::shared_ptr<BcRotVec> EshelbianPlasticity::EshelbianCore::bcSpatialRotationVecPtr
Examples
EshelbianPlasticity.cpp.

Definition at line 1216 of file EshelbianPlasticity.hpp.

◆ bcSpatialTraction

boost::shared_ptr<TractionBcVec> EshelbianPlasticity::EshelbianCore::bcSpatialTraction
Examples
EshelbianPlasticity.cpp.

Definition at line 1217 of file EshelbianPlasticity.hpp.

◆ bubbleField

const std::string EshelbianPlasticity::EshelbianCore::bubbleField
Examples
EshelbianPlasticity.cpp.

Definition at line 1197 of file EshelbianPlasticity.hpp.

◆ dataAtPts

boost::shared_ptr<DataAtIntegrationPts> EshelbianPlasticity::EshelbianCore::dataAtPts
Examples
EshelbianPlasticity.cpp.

Definition at line 1170 of file EshelbianPlasticity.hpp.

◆ dM

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dM

Coupled problem all fields.

Examples
EshelbianPlasticity.cpp.

Definition at line 1179 of file EshelbianPlasticity.hpp.

◆ dmElastic

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dmElastic

Elastic problem.

Examples
EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 1180 of file EshelbianPlasticity.hpp.

◆ dmElasticSchurBubble

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dmElasticSchurBubble

Sub problem of dmElastic Schur.

Examples
EshelbianPlasticity.cpp.

Definition at line 1182 of file EshelbianPlasticity.hpp.

◆ dmElasticSchurOmega

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dmElasticSchurOmega

Sub problem of dmElastic Schur.

Examples
EshelbianPlasticity.cpp.

Definition at line 1185 of file EshelbianPlasticity.hpp.

◆ dmElasticSchurSpatialDisp

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dmElasticSchurSpatialDisp

Sub problem of dmElastic Schur.

Examples
EshelbianPlasticity.cpp.

Definition at line 1184 of file EshelbianPlasticity.hpp.

◆ dmElasticSchurStreach

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dmElasticSchurStreach

Sub problem of dmElastic Schur.

Examples
EshelbianPlasticity.cpp.

Definition at line 1181 of file EshelbianPlasticity.hpp.

◆ dmMaterial

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dmMaterial

Material problem.

Definition at line 1186 of file EshelbianPlasticity.hpp.

◆ elasticBcLhs

boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore> > EshelbianPlasticity::EshelbianCore::elasticBcLhs
Examples
EshelbianPlasticity.cpp.

Definition at line 1175 of file EshelbianPlasticity.hpp.

◆ elasticBcRhs

boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore> > EshelbianPlasticity::EshelbianCore::elasticBcRhs
Examples
EshelbianPlasticity.cpp.

Definition at line 1176 of file EshelbianPlasticity.hpp.

◆ elasticFeLhs

boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore> > EshelbianPlasticity::EshelbianCore::elasticFeLhs
Examples
EshelbianPlasticity.cpp.

Definition at line 1174 of file EshelbianPlasticity.hpp.

◆ elasticFeRhs

boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore> > EshelbianPlasticity::EshelbianCore::elasticFeRhs
Examples
EshelbianPlasticity.cpp.

Definition at line 1173 of file EshelbianPlasticity.hpp.

◆ elementVolumeName

const std::string EshelbianPlasticity::EshelbianCore::elementVolumeName
Examples
EshelbianPlasticity.cpp.

Definition at line 1199 of file EshelbianPlasticity.hpp.

◆ eshelbyStress

const std::string EshelbianPlasticity::EshelbianCore::eshelbyStress
Examples
EshelbianPlasticity.cpp.

Definition at line 1189 of file EshelbianPlasticity.hpp.

◆ essentialBcElement

const std::string EshelbianPlasticity::EshelbianCore::essentialBcElement
Examples
EshelbianPlasticity.cpp.

Definition at line 1201 of file EshelbianPlasticity.hpp.

◆ lambdaField

const std::string EshelbianPlasticity::EshelbianCore::lambdaField

Definition at line 1196 of file EshelbianPlasticity.hpp.

◆ materialDisp

const std::string EshelbianPlasticity::EshelbianCore::materialDisp

Definition at line 1191 of file EshelbianPlasticity.hpp.

◆ materialGradient

const std::string EshelbianPlasticity::EshelbianCore::materialGradient
Examples
EshelbianPlasticity.cpp.

Definition at line 1194 of file EshelbianPlasticity.hpp.

◆ materialOrder

int EshelbianPlasticity::EshelbianCore::materialOrder
Examples
EshelbianPlasticity.cpp.

Definition at line 1207 of file EshelbianPlasticity.hpp.

◆ mField

MoFEM::Interface& EshelbianPlasticity::EshelbianCore::mField
Examples
EshelbianPlasticity.cpp.

Definition at line 1168 of file EshelbianPlasticity.hpp.

◆ naturalBcElement

const std::string EshelbianPlasticity::EshelbianCore::naturalBcElement
Examples
EshelbianPlasticity.cpp.

Definition at line 1200 of file EshelbianPlasticity.hpp.

◆ physicalEquations

boost::shared_ptr<PhysicalEquations> EshelbianPlasticity::EshelbianCore::physicalEquations
Examples
EshelbianPlasticity.cpp.

Definition at line 1171 of file EshelbianPlasticity.hpp.

◆ piolaStress

const std::string EshelbianPlasticity::EshelbianCore::piolaStress
Examples
EshelbianPlasticity.cpp.

Definition at line 1188 of file EshelbianPlasticity.hpp.

◆ precEps

double EshelbianPlasticity::EshelbianCore::precEps
Examples
EshelbianPlasticity.cpp.

Definition at line 1211 of file EshelbianPlasticity.hpp.

◆ rotAxis

const std::string EshelbianPlasticity::EshelbianCore::rotAxis
Examples
EshelbianPlasticity.cpp.

Definition at line 1193 of file EshelbianPlasticity.hpp.

◆ schurAssembly

boost::shared_ptr<EpFEMethod> EshelbianPlasticity::EshelbianCore::schurAssembly
Examples
EshelbianPlasticity.cpp.

Definition at line 1177 of file EshelbianPlasticity.hpp.

◆ spaceOrder

int EshelbianPlasticity::EshelbianCore::spaceOrder
Examples
EshelbianPlasticity.cpp.

Definition at line 1206 of file EshelbianPlasticity.hpp.

◆ spatialDisp

const std::string EshelbianPlasticity::EshelbianCore::spatialDisp
Examples
EshelbianPlasticity.cpp.

Definition at line 1190 of file EshelbianPlasticity.hpp.

◆ streachTensor

const std::string EshelbianPlasticity::EshelbianCore::streachTensor
Examples
EshelbianPlasticity.cpp.

Definition at line 1192 of file EshelbianPlasticity.hpp.

◆ tauField

const std::string EshelbianPlasticity::EshelbianCore::tauField

Definition at line 1195 of file EshelbianPlasticity.hpp.


The documentation for this struct was generated from the following files: