v0.13.2
Loading...
Searching...
No Matches
Classes | Public Member Functions | Public Attributes | List of all members
MagneticElement Struct Reference

Implementation of magneto-static problem (basic Implementation) More...

#include <users_modules/tutorials/max-0/src/MagneticElement.hpp>

Collaboration diagram for MagneticElement:
[legend]

Classes

struct  BlockData
 data structure storing material constants, model parameters, matrices, etc. More...
 
struct  OpCurlCurl
 calculate and assemble CurlCurl matrix More...
 
struct  OpNaturalBC
 calculate essential boundary conditions More...
 
struct  OpPostProcessCurl
 calculate and assemble CurlCurl matrix More...
 
struct  OpStab
 calculate and assemble stabilization matrix More...
 
struct  TriFE
 define surface element More...
 
struct  VolumeFE
 definition of volume element More...
 

Public Member Functions

 MagneticElement (MoFEM::Interface &m_field)
 
virtual ~MagneticElement ()=default
 
MoFEMErrorCode getNaturalBc ()
 get natural boundary conditions More...
 
MoFEMErrorCode getEssentialBc ()
 get essential boundary conditions (only homogenous case is considered) More...
 
MoFEMErrorCode createFields ()
 build problem data structures More...
 
MoFEMErrorCode createElements ()
 create finite elements More...
 
MoFEMErrorCode createProblem ()
 create problem More...
 
MoFEMErrorCode destroyProblem ()
 destroy Distributed mesh manager More...
 
MoFEMErrorCode solveProblem (const bool regression_test=false)
 solve problem More...
 
MoFEMErrorCode postProcessResults ()
 post-process results, i.e. save solution on the mesh More...
 

Public Attributes

MoFEM::InterfacemField
 
BlockData blockData
 

Detailed Description

Implementation of magneto-static problem (basic Implementation)

Look for theory and details here:

[30] <www.hpfem.jku.at/publications/szthesis.pdf> https://pdfs.semanticscholar.org/0574/e5763d9b64bff16f908e9621f23af3b3dc86.pdf

Election file and all other problem related file are here maxwell_element.

Todo:

Extension for mix formulation

Use appropriate pre-conditioner for large problems

Examples
MagneticElement.hpp, and magnetostatic.cpp.

Definition at line 28 of file MagneticElement.hpp.

Constructor & Destructor Documentation

◆ MagneticElement()

MagneticElement::MagneticElement ( MoFEM::Interface m_field)
inline

Definition at line 56 of file MagneticElement.hpp.

56: mField(m_field) {}
MoFEM::Interface & mField

◆ ~MagneticElement()

virtual MagneticElement::~MagneticElement ( )
virtualdefault

Member Function Documentation

◆ createElements()

MoFEMErrorCode MagneticElement::createElements ( )
inline

create finite elements

Create volume and surface element. Surface element is used to integrate natural boundary conditions.

Returns
error code
Examples
MagneticElement.hpp, and magnetostatic.cpp.

Definition at line 234 of file MagneticElement.hpp.

234 {
236 // //Elements
245 "MESH_NODE_POSITIONS");
254 blockData.feNaturalBCName, "MESH_NODE_POSITIONS");
259 // build finite elemnts
261 // build adjacencies
264 }
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#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
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies

◆ createFields()

MoFEMErrorCode MagneticElement::createFields ( )
inline

build problem data structures

Returns
error code
Examples
MagneticElement.hpp, and magnetostatic.cpp.

Definition at line 158 of file MagneticElement.hpp.

158 {
160
161 // Set entities bit level. each entity has bit level depending for example
162 // on refinement level. In this case we do not refine mesh or not do
163 // topological changes, simply set refinement level to zero on all entities.
164
165 CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(
166 0, 3, BitRefLevel().set(0));
167
168 // add fields
170 1);
171 CHKERR mField.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
172 3);
173 // meshset consisting all entities in mesh
174 EntityHandle root_set = mField.get_moab().get_root_set();
175 // add entities to field
178
179 // // The higher-order gradients can be gauged by locally skipping the
180 // // corresponding degrees of freedom and basis functions in the
181 // higher-order
182 // // edge-based, face-based and cell-based finite element subspaces.
183 //
184 // Range tris,edges;
185 // CHKERR mField.get_moab().get_entities_by_type(root_set,MBTRI,tris,true);
186 // mField.get_moab().get_entities_by_type(root_set,MBEDGE,edges,true);
187 //
188 // // Set order in volume
189 // Range bc_ents = unite(blockData.naturalBc,blockData.essentialBc);
190 // Range vol_ents = subtract(unite(tris,edges),bc_ents);
191 // CHKERR
192 // mField.set_field_order(vol_ents,blockData.fieldName,blockData.oRder); int
193 // gauged_order = 1; CHKERR
194 // mField.set_field_order(bc_ents,blockData.fieldName,gauged_order);
195
196 // Set order on tets
203
204 // Set geometry approximation ordered
206 "MESH_NODE_POSITIONS");
207 CHKERR mField.set_field_order(root_set, MBTET, "MESH_NODE_POSITIONS", 2);
208 CHKERR mField.set_field_order(root_set, MBTRI, "MESH_NODE_POSITIONS", 2);
209 CHKERR mField.set_field_order(root_set, MBEDGE, "MESH_NODE_POSITIONS", 2);
210 CHKERR mField.set_field_order(root_set, MBVERTEX, "MESH_NODE_POSITIONS", 1);
211
212 // build field
214
215 // get HO geometry for 10 node tets
216 // This method takes coordinates form edges mid nodes in 10 node tet and
217 // project values on 2nd order hierarchical basis used to approx. geometry.
218 Projection10NodeCoordsOnField ent_method_material(mField,
219 "MESH_NODE_POSITIONS");
220 CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
221
223 }
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
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 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.
int oRder
approximation order
virtual moab::Interface & get_moab()=0
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.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

◆ createProblem()

MoFEMErrorCode MagneticElement::createProblem ( )
inline

create problem

Problem is collection of finite elements. With the information on which fields finite elements operates the matrix and left and right hand side vector could be created.

Here we use Distributed mesh manager from PETSc as a inteface.

Returns
error code
Examples
MagneticElement.hpp, and magnetostatic.cpp.

Definition at line 278 of file MagneticElement.hpp.

278 {
280 // set up DM
281 CHKERR DMRegister_MoFEM("DMMOFEM");
282 CHKERR DMCreate(PETSC_COMM_WORLD, &blockData.dM);
283 CHKERR DMSetType(blockData.dM, "DMMOFEM");
284 CHKERR DMMoFEMCreateMoFEM(blockData.dM, &mField, "MAGNETIC_PROBLEM",
285 BitRefLevel().set(0));
286 CHKERR DMSetFromOptions(blockData.dM);
288 // add elements to blockData.dM
291 CHKERR DMSetUp(blockData.dM);
292
293 // remove essential DOFs
294 const MoFEM::Problem *problem_ptr;
296 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
298
300 }
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1091
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:465
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 DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:394
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
keeps basic data about problem

◆ destroyProblem()

MoFEMErrorCode MagneticElement::destroyProblem ( )
inline

destroy Distributed mesh manager

Returns
[description]
Examples
MagneticElement.hpp, and magnetostatic.cpp.

Definition at line 307 of file MagneticElement.hpp.

307 {
309 CHKERR DMDestroy(&blockData.dM);
311 }

◆ getEssentialBc()

MoFEMErrorCode MagneticElement::getEssentialBc ( )
inline

get essential boundary conditions (only homogenous case is considered)

Returns
error code
Examples
MagneticElement.hpp, and magnetostatic.cpp.

Definition at line 121 of file MagneticElement.hpp.

121 {
122 ParallelComm *pcomm =
123 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
126 if (bit->getName().compare(0, 10, "ESSENTIALBC") == 0) {
127 Range faces;
128 CHKERR mField.get_moab().get_entities_by_type(bit->meshset, MBTRI,
129 faces, true);
130 CHKERR mField.get_moab().get_adjacencies(
131 faces, 1, true, blockData.essentialBc, moab::Interface::UNION);
132 blockData.essentialBc.merge(faces);
133 }
134 }
135 if (blockData.essentialBc.empty()) {
136 Range tets;
137 CHKERR mField.get_moab().get_entities_by_type(0, MBTET, tets);
138 Skinner skin(&mField.get_moab());
139 Range skin_faces; // skin faces from 3d ents
140 CHKERR skin.find_skin(0, tets, false, skin_faces);
141 skin_faces = subtract(skin_faces, blockData.naturalBc);
142 Range proc_skin;
143 CHKERR pcomm->filter_pstatus(skin_faces,
144 PSTATUS_SHARED | PSTATUS_MULTISHARED,
145 PSTATUS_NOT, -1, &proc_skin);
146 CHKERR mField.get_moab().get_adjacencies(
147 proc_skin, 1, true, blockData.essentialBc, moab::Interface::UNION);
148 blockData.essentialBc.merge(proc_skin);
149 }
151 }
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
@ BLOCKSET
Definition: definitions.h:148
#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.
auto bit
set bit

◆ getNaturalBc()

MoFEMErrorCode MagneticElement::getNaturalBc ( )
inline

get natural boundary conditions

Returns
error code
Examples
MagneticElement.hpp, and magnetostatic.cpp.

Definition at line 101 of file MagneticElement.hpp.

101 {
104 if (bit->getName().compare(0, 9, "NATURALBC") == 0) {
105 Range faces;
106 CHKERR mField.get_moab().get_entities_by_type(bit->meshset, MBTRI,
107 faces, true);
108 CHKERR mField.get_moab().get_adjacencies(
109 faces, 1, true, blockData.naturalBc, moab::Interface::UNION);
110 blockData.naturalBc.merge(faces);
111 }
112 }
114 }

◆ postProcessResults()

MoFEMErrorCode MagneticElement::postProcessResults ( )
inline

post-process results, i.e. save solution on the mesh

Returns
[description]
Examples
MagneticElement.hpp, and magnetostatic.cpp.

Definition at line 406 of file MagneticElement.hpp.

406 {
407
409 PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(
410 mField);
411
412 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", post_proc, false, true, false,
413 true);
414
415 auto pos_ptr = boost::make_shared<MatrixDouble>();
416 auto field_val_ptr = boost::make_shared<MatrixDouble>();
417
418 post_proc.getOpPtrVector().push_back(
419 new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS", pos_ptr));
420 post_proc.getOpPtrVector().push_back(
421 new OpCalculateHVecVectorField<3>(blockData.fieldName, field_val_ptr));
422
423 using OpPPMap = OpPostProcMapInMoab<3, 3>;
424
425 post_proc.getOpPtrVector().push_back(
426
427 new OpPPMap(
428
429 post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
430
431 {},
432
433 {{"MESH_NODE_POSITIONS", pos_ptr},
434 {blockData.fieldName, field_val_ptr}},
435
436 {},
437
438 {}
439
440 )
441
442 );
443
444 post_proc.getOpPtrVector().push_back(new OpPostProcessCurl(
445 blockData, post_proc.getPostProcMesh(), post_proc.getMapGaussPts()));
447 &post_proc);
448 CHKERR post_proc.writeFile("out_values.h5m");
450 }
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:554
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Post post-proc data at points from hash maps.

◆ solveProblem()

MoFEMErrorCode MagneticElement::solveProblem ( const bool  regression_test = false)
inline

solve problem

Create matrices; integrate over elements; solve linear system of equations

Examples
MagneticElement.hpp, and magnetostatic.cpp.

Definition at line 319 of file MagneticElement.hpp.

319 {
321
322 VolumeFE vol_fe(mField);
323 auto material_grad_mat = boost::make_shared<MatrixDouble>();
324 auto material_det_vec = boost::make_shared<VectorDouble>();
325 auto material_inv_grad_mat = boost::make_shared<MatrixDouble>();
326 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", vol_fe, false, true, false, true);
327 vol_fe.getOpPtrVector().push_back(new OpCurlCurl(blockData));
328 vol_fe.getOpPtrVector().push_back(new OpStab(blockData));
329 TriFE tri_fe(mField);
330 tri_fe.meshPositionsFieldName = "none";
331
332 tri_fe.getOpPtrVector().push_back(
333 new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
334 tri_fe.getOpPtrVector().push_back(
335 new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
336 tri_fe.getOpPtrVector().push_back(
337 new OpHOSetCovariantPiolaTransformOnFace3D(HCURL));
338 tri_fe.getOpPtrVector().push_back(new OpNaturalBC(blockData));
339
340 // create matrices and vectors
341 CHKERR DMCreateGlobalVector(blockData.dM, &blockData.D);
342 // CHKERR DMCreateMatrix(blockData.dM, &blockData.A);
343 CHKERR DMCreateGlobalVector(blockData.dM, &blockData.F);
344 CHKERR VecDuplicate(blockData.F, &blockData.D);
345
346 const MoFEM::Problem *problem_ptr;
348 CHKERR mField.getInterface<MatrixManager>()
349 ->createMPIAIJ<PetscGlobalIdx_mi_tag>(problem_ptr->getName(),
350 &blockData.A);
351
352 CHKERR MatZeroEntries(blockData.A);
353 CHKERR VecZeroEntries(blockData.F);
354 CHKERR VecGhostUpdateBegin(blockData.F, ADD_VALUES, SCATTER_FORWARD);
355 CHKERR VecGhostUpdateEnd(blockData.F, ADD_VALUES, SCATTER_FORWARD);
356
358 &vol_fe);
360 blockData.feNaturalBCName.c_str(), &tri_fe);
361
362 CHKERR MatAssemblyBegin(blockData.A, MAT_FINAL_ASSEMBLY);
363 CHKERR MatAssemblyEnd(blockData.A, MAT_FINAL_ASSEMBLY);
364 CHKERR VecGhostUpdateBegin(blockData.F, ADD_VALUES, SCATTER_REVERSE);
365 CHKERR VecGhostUpdateEnd(blockData.F, ADD_VALUES, SCATTER_REVERSE);
366 CHKERR VecAssemblyBegin(blockData.F);
367 CHKERR VecAssemblyEnd(blockData.F);
368
369 KSP solver;
370 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
371 CHKERR KSPSetOperators(solver, blockData.A, blockData.A);
372 CHKERR KSPSetFromOptions(solver);
373 CHKERR KSPSetUp(solver);
374 CHKERR KSPSolve(solver, blockData.F, blockData.D);
375 CHKERR KSPDestroy(&solver);
376
377 CHKERR VecGhostUpdateBegin(blockData.D, INSERT_VALUES, SCATTER_FORWARD);
378 CHKERR VecGhostUpdateEnd(blockData.D, INSERT_VALUES, SCATTER_FORWARD);
380 SCATTER_REVERSE);
381
382 if (regression_test) {
383 // This test is for order = 1 only
384 double nrm2;
385 CHKERR VecNorm(blockData.D, NORM_2, &nrm2);
386 const double expected_value = 4.6772e+01;
387 const double tol = 1e-4;
388 if ((nrm2 - expected_value) / expected_value > tol) {
389 SETERRQ2(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
390 "Regression test failed %6.4e != %6.4e", nrm2, expected_value);
391 }
392 }
393
394 CHKERR VecDestroy(&blockData.D);
395 CHKERR VecDestroy(&blockData.F);
396 CHKERR MatDestroy(&blockData.A);
397
399 }
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
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:491
double tol

Member Data Documentation

◆ blockData

BlockData MagneticElement::blockData
Examples
MagneticElement.hpp, and magnetostatic.cpp.

Definition at line 94 of file MagneticElement.hpp.

◆ mField

MoFEM::Interface& MagneticElement::mField
Examples
MagneticElement.hpp.

Definition at line 30 of file MagneticElement.hpp.


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