MagneticElement Struct Reference

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

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

Collaboration diagram for MagneticElement:

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::Interface &	mField
BlockData	blockData

Detailed Description

Implementation of magneto-static problem (basic Implementation)

Look for theory and details here:

[35] <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

Examples

MagneticElement.hpp.

Definition at line 56 of file MagneticElement.hpp.

: mField(m_field) {}

~MagneticElement()

virtual MagneticElement::~MagneticElement ( )

virtualdefault

Examples

MagneticElement.hpp.

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.

                                  {
    MoFEMFunctionBegin;
    // //Elements
    CHKERR mField.add_finite_element(blockData.feName);
    CHKERR mField.modify_finite_element_add_field_row(blockData.feName,
                                                      blockData.fieldName);
    CHKERR mField.modify_finite_element_add_field_col(blockData.feName,
                                                      blockData.fieldName);
    CHKERR mField.modify_finite_element_add_field_data(blockData.feName,
                                                       blockData.fieldName);
    CHKERR mField.modify_finite_element_add_field_data(blockData.feName,
                                                       "MESH_NODE_POSITIONS");
    CHKERR mField.add_finite_element(blockData.feNaturalBCName);
    CHKERR mField.modify_finite_element_add_field_row(blockData.feNaturalBCName,
                                                      blockData.fieldName);
    CHKERR mField.modify_finite_element_add_field_col(blockData.feNaturalBCName,
                                                      blockData.fieldName);
    CHKERR mField.modify_finite_element_add_field_data(
        blockData.feNaturalBCName, blockData.fieldName);
    CHKERR mField.modify_finite_element_add_field_data(
        blockData.feNaturalBCName, "MESH_NODE_POSITIONS");
    CHKERR mField.add_ents_to_finite_element_by_type(0, MBTET,
                                                     blockData.feName);
    CHKERR mField.add_ents_to_finite_element_by_type(blockData.naturalBc, MBTRI,
                                                     blockData.feNaturalBCName);
    // build finite elemnts
    CHKERR mField.build_finite_elements();
    // build adjacencies
    CHKERR mField.build_adjacencies(BitRefLevel().set(0));
    MoFEMFunctionReturn(0);
  }

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.

                                {
    MoFEMFunctionBegin;
 
    // Set entities bit level. each entity has bit level depending for example
    // on refinement level. In this case we do not refine mesh or not do
    // topological changes, simply set refinement level to zero on all entities.
 
    CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(
        0, 3, BitRefLevel().set(0));
 
    // add fields
    CHKERR mField.add_field(blockData.fieldName, HCURL, DEMKOWICZ_JACOBI_BASE,
                            1);
    CHKERR mField.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
                            3);
    // meshset consisting all entities in mesh
    EntityHandle root_set = mField.get_moab().get_root_set();
    // add entities to field
    CHKERR mField.add_ents_to_field_by_type(root_set, MBTET,
                                            blockData.fieldName);
 
    // // The higher-order gradients can be gauged by locally skipping the
    // // corresponding degrees of freedom and basis functions in the
    // higher-order
    // // edge-based, face-based and cell-based finite element subspaces.
    //
    // Range tris,edges;
    // CHKERR mField.get_moab().get_entities_by_type(root_set,MBTRI,tris,true);
    // mField.get_moab().get_entities_by_type(root_set,MBEDGE,edges,true);
    //
    // // Set order in volume
    // Range bc_ents = unite(blockData.naturalBc,blockData.essentialBc);
    // Range vol_ents = subtract(unite(tris,edges),bc_ents);
    // CHKERR
    // mField.set_field_order(vol_ents,blockData.fieldName,blockData.oRder); int
    // gauged_order = 1; CHKERR
    // mField.set_field_order(bc_ents,blockData.fieldName,gauged_order);
 
    // Set order on tets
    CHKERR mField.set_field_order(root_set, MBTET, blockData.fieldName,
                                  blockData.oRder);
    CHKERR mField.set_field_order(root_set, MBTRI, blockData.fieldName,
                                  blockData.oRder);
    CHKERR mField.set_field_order(root_set, MBEDGE, blockData.fieldName,
                                  blockData.oRder);
 
    // Set geometry approximation ordered
    CHKERR mField.add_ents_to_field_by_type(root_set, MBTET,
                                            "MESH_NODE_POSITIONS");
    CHKERR mField.set_field_order(root_set, MBTET, "MESH_NODE_POSITIONS", 2);
    CHKERR mField.set_field_order(root_set, MBTRI, "MESH_NODE_POSITIONS", 2);
    CHKERR mField.set_field_order(root_set, MBEDGE, "MESH_NODE_POSITIONS", 2);
    CHKERR mField.set_field_order(root_set, MBVERTEX, "MESH_NODE_POSITIONS", 1);
 
    // build field
    CHKERR mField.build_fields();
 
    // get HO geometry for 10 node tets
    // This method takes coordinates form edges mid nodes in 10 node tet and
    // project values on 2nd order hierarchical basis used to approx. geometry.
    Projection10NodeCoordsOnField ent_method_material(mField,
                                                      "MESH_NODE_POSITIONS");
    CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
 
    MoFEMFunctionReturn(0);
  }

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.

                                 {
    MoFEMFunctionBegin;
    // set up DM
    CHKERR DMRegister_MoFEM("DMMOFEM");
    CHKERR DMCreate(PETSC_COMM_WORLD, &blockData.dM);
    CHKERR DMSetType(blockData.dM, "DMMOFEM");
    CHKERR DMMoFEMCreateMoFEM(blockData.dM, &mField, "MAGNETIC_PROBLEM",
                              BitRefLevel().set(0));
    CHKERR DMSetFromOptions(blockData.dM);
    CHKERR DMMoFEMSetIsPartitioned(blockData.dM, PETSC_TRUE);
    // add elements to blockData.dM
    CHKERR DMMoFEMAddElement(blockData.dM, blockData.feName);
    CHKERR DMMoFEMAddElement(blockData.dM, blockData.feNaturalBCName);
    CHKERR DMSetUp(blockData.dM);
 
    // remove essential DOFs
    const MoFEM::Problem *problem_ptr;
    CHKERR DMMoFEMGetProblemPtr(blockData.dM, &problem_ptr);
    CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
        problem_ptr->getName(), blockData.fieldName, blockData.essentialBc);
 
    MoFEMFunctionReturn(0);
  }

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.

                                  {
    MoFEMFunctionBegin;
    CHKERR DMDestroy(&blockData.dM);
    MoFEMFunctionReturn(0);
  }

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.

                                  {
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
    MoFEMFunctionBegin;
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, bit)) {
      if (bit->getName().compare(0, 10, "ESSENTIALBC") == 0) {
        Range faces;
        CHKERR mField.get_moab().get_entities_by_type(bit->meshset, MBTRI,
                                                      faces, true);
        CHKERR mField.get_moab().get_adjacencies(
            faces, 1, true, blockData.essentialBc, moab::Interface::UNION);
        blockData.essentialBc.merge(faces);
      }
    }
    if (blockData.essentialBc.empty()) {
      Range tets;
      CHKERR mField.get_moab().get_entities_by_type(0, MBTET, tets);
      Skinner skin(&mField.get_moab());
      Range skin_faces; // skin faces from 3d ents
      CHKERR skin.find_skin(0, tets, false, skin_faces);
      skin_faces = subtract(skin_faces, blockData.naturalBc);
      Range proc_skin;
      CHKERR pcomm->filter_pstatus(skin_faces,
                                   PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                   PSTATUS_NOT, -1, &proc_skin);
      CHKERR mField.get_moab().get_adjacencies(
          proc_skin, 1, true, blockData.essentialBc, moab::Interface::UNION);
      blockData.essentialBc.merge(proc_skin);
    }
    MoFEMFunctionReturn(0);
  }

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.

                                {
    MoFEMFunctionBegin;
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, bit)) {
      if (bit->getName().compare(0, 9, "NATURALBC") == 0) {
        Range faces;
        CHKERR mField.get_moab().get_entities_by_type(bit->meshset, MBTRI,
                                                      faces, true);
        CHKERR mField.get_moab().get_adjacencies(
            faces, 1, true, blockData.naturalBc, moab::Interface::UNION);
        blockData.naturalBc.merge(faces);
      }
    }
    MoFEMFunctionReturn(0);
  }

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.

                                      {
 
    MoFEMFunctionBegin;
    PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(
        mField);
 
    CHKERR addHOOpsVol("MESH_NODE_POSITIONS", post_proc, false, true, false,
                       true);
 
    auto pos_ptr = boost::make_shared<MatrixDouble>();
    auto field_val_ptr = boost::make_shared<MatrixDouble>();
 
    post_proc.getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS", pos_ptr));
    post_proc.getOpPtrVector().push_back(
        new OpCalculateHVecVectorField<3>(blockData.fieldName, field_val_ptr));
 
    using OpPPMap = OpPostProcMapInMoab<3, 3>;
 
    post_proc.getOpPtrVector().push_back(
 
        new OpPPMap(
 
            post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
 
            {},
 
            {{"MESH_NODE_POSITIONS", pos_ptr},
             {blockData.fieldName, field_val_ptr}},
 
            {},
 
            {}
 
            )
 
    );
 
    post_proc.getOpPtrVector().push_back(new OpPostProcessCurl(
        blockData, post_proc.getPostProcMesh(), post_proc.getMapGaussPts()));
    CHKERR DMoFEMLoopFiniteElements(blockData.dM, blockData.feName.c_str(),
                                    &post_proc);
    CHKERR post_proc.writeFile("out_values.h5m");
    MoFEMFunctionReturn(0);
  }

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.

                                                                  {
    MoFEMFunctionBegin;
 
    VolumeFE vol_fe(mField);
    auto material_grad_mat = boost::make_shared<MatrixDouble>();
    auto material_det_vec = boost::make_shared<VectorDouble>();
    auto material_inv_grad_mat = boost::make_shared<MatrixDouble>();
    CHKERR addHOOpsVol("MESH_NODE_POSITIONS", vol_fe, false, true, false, true);
    vol_fe.getOpPtrVector().push_back(new OpCurlCurl(blockData));
    vol_fe.getOpPtrVector().push_back(new OpStab(blockData));
    TriFE tri_fe(mField);
    tri_fe.meshPositionsFieldName = "none";
 
    tri_fe.getOpPtrVector().push_back(
        new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
    tri_fe.getOpPtrVector().push_back(
        new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
    tri_fe.getOpPtrVector().push_back(
        new OpHOSetCovariantPiolaTransformOnFace3D(HCURL));
    tri_fe.getOpPtrVector().push_back(new OpNaturalBC(blockData));
 
    // create matrices and vectors
    CHKERR DMCreateGlobalVector(blockData.dM, &blockData.D);
    // CHKERR DMCreateMatrix(blockData.dM, &blockData.A);
    CHKERR DMCreateGlobalVector(blockData.dM, &blockData.F);
    CHKERR VecDuplicate(blockData.F, &blockData.D);
 
    const MoFEM::Problem *problem_ptr;
    CHKERR DMMoFEMGetProblemPtr(blockData.dM, &problem_ptr);
    CHKERR mField.getInterface<MatrixManager>()
        ->createMPIAIJ<PetscGlobalIdx_mi_tag>(problem_ptr->getName(),
                                              &blockData.A);
 
    CHKERR MatZeroEntries(blockData.A);
    CHKERR VecZeroEntries(blockData.F);
    CHKERR VecGhostUpdateBegin(blockData.F, ADD_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(blockData.F, ADD_VALUES, SCATTER_FORWARD);
 
    CHKERR DMoFEMLoopFiniteElements(blockData.dM, blockData.feName.c_str(),
                                    &vol_fe);
    CHKERR DMoFEMLoopFiniteElements(blockData.dM,
                                    blockData.feNaturalBCName.c_str(), &tri_fe);
 
    CHKERR MatAssemblyBegin(blockData.A, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(blockData.A, MAT_FINAL_ASSEMBLY);
    CHKERR VecGhostUpdateBegin(blockData.F, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(blockData.F, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecAssemblyBegin(blockData.F);
    CHKERR VecAssemblyEnd(blockData.F);
 
    KSP solver;
    CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
    CHKERR KSPSetOperators(solver, blockData.A, blockData.A);
    CHKERR KSPSetFromOptions(solver);
    CHKERR KSPSetUp(solver);
    CHKERR KSPSolve(solver, blockData.F, blockData.D);
    CHKERR KSPDestroy(&solver);
 
    CHKERR VecGhostUpdateBegin(blockData.D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(blockData.D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR DMoFEMMeshToLocalVector(blockData.dM, blockData.D, INSERT_VALUES,
                                   SCATTER_REVERSE);
 
    if (regression_test) {
      // This test is for order = 1 only
      double nrm2;
      CHKERR VecNorm(blockData.D, NORM_2, &nrm2);
      const double expected_value = 4.6772e+01;
      const double tol = 1e-4;
      if ((nrm2 - expected_value) / expected_value > tol) {
        SETERRQ2(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                 "Regression test failed %6.4e != %6.4e", nrm2, expected_value);
      }
    }
 
    CHKERR VecDestroy(&blockData.D);
    CHKERR VecDestroy(&blockData.F);
    CHKERR MatDestroy(&blockData.A);
 
    MoFEMFunctionReturn(0);
  }

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:

• tutorials/max-0/src/MagneticElement.hpp