MagneticElement Struct Reference

Implementation of magneto-static problem (basic Implementation)Look for theory and details here: More...

#include <users_modules/basic_finite_elements/magnetostatic/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 \[ \mathbf{A} = \int_\Omega \mu^{-1} \left( \nabla \times \mathbf{u} \cdot \nabla \times \mathbf{v} \right) \textrm{d}\Omega \] where \[ \mathbf{u} = \nabla \times \mathbf{B} \] where \(\mathbf{B}\) is magnetic flux and \(\mu\) is magnetic permeability. More...
struct	OpNaturalBC
	calculate essential boundary conditions \[ \mathbf{F} = \int_{\partial\Omega} \ \mathbf{u} \cdot \mathbf{j}_i \textrm{d}{\partial\Omega} \] where \(\mathbf{j}_i\) is current density function. More...
struct	OpPostProcessCurl
	calculate and assemble CurlCurl matrix More...
struct	OpStab
	calculate and assemble stabilization matrix \[ \mathbf{A} = \int_\Omega \epsilon \mathbf{u} \cdot \mathbf{v} \textrm{d}\Omega \] where \(\epsilon\) is regularization parameter. More...
struct	TriFE
	define surface element More...
struct	VolumeFE
	definition of volume element More...

Public Member Functions
	MagneticElement (MoFEM::Interface &m_field)
virtual	~MagneticElement ()
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 elementsCreate volume and surface element. Surface element is used to integrate natural boundary conditions. More...
MoFEMErrorCode	createProblem ()
	create problemProblem 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. More...
MoFEMErrorCode	destroyProblem ()
	destroy Distributed mesh manager More...
MoFEMErrorCode	solveProblem (const bool regression_test=false)
	solve problemCreate matrices; integrate over elements; solve linear system of equations 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:

[28] <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 43 of file MagneticElement.hpp.

Constructor & Destructor Documentation

MagneticElement()

MagneticElement::MagneticElement

(

MoFEM::Interface &

m_field

)

Examples

MagneticElement.hpp.

Definition at line 71 of file MagneticElement.hpp.

: mField(m_field) {}

~MagneticElement()

virtual MagneticElement::~MagneticElement ( )

virtual

Examples

MagneticElement.hpp.

Definition at line 72 of file MagneticElement.hpp.

{}

Member Function Documentation

createElements()

MoFEMErrorCode MagneticElement::createElements

(

)

create finite elementsCreate 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 243 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

(

)

build problem data structures

Returns

error code

Examples

MagneticElement.hpp, and magnetostatic.cpp.

Definition at line 167 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

(

)

create problemProblem 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 287 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);
    // add elements to blockData.dM
    CHKERR DMMoFEMAddElement(blockData.dM, blockData.feName.c_str());
    CHKERR DMMoFEMAddElement(blockData.dM, blockData.feNaturalBCName.c_str());
    CHKERR DMSetUp(blockData.dM);
    // create matrices and vectors
    CHKERR DMCreateGlobalVector(blockData.dM, &blockData.D);
    MoFEMFunctionReturn(0);
  }

destroyProblem()

MoFEMErrorCode MagneticElement::destroyProblem

(

)

destroy Distributed mesh manager

Returns

[description]

Examples

MagneticElement.hpp, and magnetostatic.cpp.

Definition at line 310 of file MagneticElement.hpp.

                                  {

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

getEssentialBc()

MoFEMErrorCode MagneticElement::getEssentialBc

(

)

get essential boundary conditions (only homogenous case is considered)

Returns

error code

Examples

MagneticElement.hpp, and magnetostatic.cpp.

Definition at line 136 of file MagneticElement.hpp.

                                  {
    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);
      CHKERR mField.get_moab().get_adjacencies(
          skin_faces, 1, true, blockData.essentialBc, moab::Interface::UNION);
      blockData.essentialBc.merge(skin_faces);
    }
    MoFEMFunctionReturn(0);
  }

getNaturalBc()

MoFEMErrorCode MagneticElement::getNaturalBc

(

)

get natural boundary conditions

Returns

error code

Examples

MagneticElement.hpp, and magnetostatic.cpp.

Definition at line 116 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

(

)

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

Returns

[description]

Examples

MagneticElement.hpp, and magnetostatic.cpp.

Definition at line 412 of file MagneticElement.hpp.

                                      {

    MoFEMFunctionBegin;
    PostProcVolumeOnRefinedMesh post_proc(mField);
    CHKERR post_proc.generateReferenceElementMesh();
    CHKERR post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
    CHKERR post_proc.addFieldValuesPostProc(blockData.fieldName);
    post_proc.getOpPtrVector().push_back(new OpPostProcessCurl(
        blockData, post_proc.postProcMesh, post_proc.mapGaussPts));
    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

)

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

Examples

MagneticElement.hpp, and magnetostatic.cpp.

Definition at line 323 of file MagneticElement.hpp.

                                                                  {
    MoFEMFunctionBegin;
    CHKERR DMCreateMatrix(blockData.dM, &blockData.A);
    CHKERR DMCreateGlobalVector(blockData.dM, &blockData.F);
    CHKERR VecDuplicate(blockData.F, &blockData.D);

    VolumeFE vol_fe(mField);
    vol_fe.getOpPtrVector().push_back(new OpCurlCurl(blockData));
    vol_fe.getOpPtrVector().push_back(new OpStab(blockData));
    TriFE tri_fe(mField);
    tri_fe.getOpPtrVector().push_back(new OpNaturalBC(blockData));

    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);

    // Boundary conditions
    std::vector<int> dofs_bc_indices;
    const MoFEM::Problem *problem_ptr;
    ierr = DMMoFEMGetProblemPtr(blockData.dM, &problem_ptr);
    for (Range::iterator eit = blockData.essentialBc.begin();
         eit != blockData.essentialBc.end(); eit++) {
      for (_IT_NUMEREDDOF_ROW_BY_NAME_ENT_PART_FOR_LOOP_(
               problem_ptr, blockData.fieldName, *eit, mField.get_comm_rank(),
               dof_ptr)) {
        std::bitset<8> pstatus(dof_ptr->get()->getPStatus());
        if (pstatus.test(0))
          continue; // only local
        dofs_bc_indices.push_back(dof_ptr->get()->getPetscGlobalDofIdx());
      }
    }

    const double diag = 1;
    CHKERR MatZeroRowsColumns(
        blockData.A, dofs_bc_indices.size(),
        dofs_bc_indices.empty() ? PETSC_NULL : &*dofs_bc_indices.begin(), diag,
        PETSC_NULL, PETSC_NULL);

    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 VecGhostUpdateBegin(blockData.D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(blockData.D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR KSPDestroy(&solver);
    CHKERR VecDestroy(&blockData.F);
    CHKERR MatDestroy(&blockData.A);

    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);
    MoFEMFunctionReturn(0);
  }

Member Data Documentation

blockData

BlockData MagneticElement::blockData

Examples

MagneticElement.hpp, and magnetostatic.cpp.

Definition at line 109 of file MagneticElement.hpp.

mField

MoFEM::Interface& MagneticElement::mField

Examples

MagneticElement.hpp.

Definition at line 45 of file MagneticElement.hpp.

---

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

• users_modules/basic_finite_elements/magnetostatic/src/MagneticElement.hpp