v0.14.0
MAX-0: Magnetostatics
Note
Prerequisites of this tutorial include SCL-1: Poisson's equation (homogeneous BC)

Note
Intended learning outcome:
• general structure of a program developed using MoFEM
• working with vectorial base fields like fields for H-curl space
• how to push the developed UDOs to the finite element
• utilisation of tools to convert outputs (MOAB) and visualise them (Paraview)
/** \file magnetostatic.cpp
* \example magnetostatic.cpp
* \ingroup maxwell_element
*
* \brief Example implementation of magnetostatics problem
*
*/
using namespace MoFEM;
static char help[] = "-my_file mesh file name\n"
"-my_order default approximation order\n"
"-my_is_partitioned set if mesh is partitioned\n"
"-regression_test check solution vector against expected value\n"
"\n\n";
int main(int argc, char *argv[]) {
const string default_options = "-ksp_type fgmres \n"
"-pc_type lu \n"
"-pc_factor_mat_solver_type mumps \n"
"-mat_mumps_icntl_20 0 \n"
"-mat_mumps_icntl_13 1 \n"
"-ksp_monitor \n"
"-mat_mumps_icntl_24 1 \n"
"-mat_mumps_icntl_13 1 \n";
string param_file = "param_file.petsc";
if (!static_cast<bool>(ifstream(param_file))) {
std::ofstream file(param_file.c_str(), std::ios::ate);
if (file.is_open()) {
file << default_options;
file.close();
}
}
MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
try {
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
auto moab_comm_wrap =
boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
if (pcomm == NULL)
pcomm =
new ParallelComm(&moab, moab_comm_wrap->get_comm());
// Read parameters from line command
PetscBool flg_file;
char mesh_file_name[255];
PetscInt order = 2;
PetscBool is_partitioned = PETSC_FALSE;
PetscBool regression_test = PETSC_FALSE;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Magnetostatics options",
"none");
CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
mesh_file_name, 255, &flg_file);
CHKERR PetscOptionsInt("-my_order", "default approximation order", "",
order, &order, PETSC_NULL);
CHKERR PetscOptionsBool(
"-regression_test",
"if set norm of solution vector is check agains expected value ",
"",
PETSC_FALSE, &regression_test, PETSC_NULL);
ierr = PetscOptionsEnd();
const char *option;
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION;";
// Create mofem interface
MoFEM::Core core(moab);
MoFEM::Interface &m_field = core;
MagneticElement magnetic(m_field);
magnetic.blockData.oRder = order;
CHKERR magnetic.getNaturalBc();
CHKERR magnetic.getEssentialBc();
CHKERR magnetic.createFields();
CHKERR magnetic.createElements();
CHKERR magnetic.createProblem();
CHKERR magnetic.solveProblem(regression_test == PETSC_TRUE);
CHKERR magnetic.postProcessResults();
CHKERR magnetic.destroyProblem();
// write solution to file (can be used by lorentz_force example)
CHKERR moab.write_file("solution.h5m", "MOAB", "PARALLEL=WRITE_PART");
}
return 0;
}
/** \file MagneticElement.hpp
* \brief Implementation of magnetic element
* \ingroup maxwell_element
* \example MagneticElement.hpp
*
*/
#ifndef __MAGNETICELEMENT_HPP__
#define __MAGNETICELEMENT_HPP__
/**
* \brief Implementation of magneto-static problem (basic Implementation)
* \ingroup maxwell_element
*
* Look for theory and details here:
*
* \cite ivanyshyn2013computation
* <www.hpfem.jku.at/publications/szthesis.pdf>
* <https://pdfs.semanticscholar.org/0574/e5763d9b64bff16f908e9621f23af3b3dc86.pdf>
*
* Election file and all other problem related file are here \ref
* maxwell_element.
*
* \todo Extension for mix formulation
* \todo Use appropriate pre-conditioner for large problems
*
*/
/// \brief definition of volume element
struct VolumeFE : public MoFEM::VolumeElementForcesAndSourcesCore {
int getRule(int order) { return 2 * order + 1; };
};
// /// \brief definition of volume element
// struct VolumeFEReducedIntegration: public
// MoFEM::VolumeElementForcesAndSourcesCore {
// VolumeFEReducedIntegration(MoFEM::Interface &m_field):
// MoFEM::VolumeElementForcesAndSourcesCore(m_field) {}
// int getRule(int order) { return 2*order+1; };
// };
/** \brief define surface element
*
*/
int getRule(int order) { return 2 * order; };
};
MagneticElement(MoFEM::Interface &m_field) : mField(m_field) {}
virtual ~MagneticElement() = default;
/**
* \brief data structure storing material constants, model parameters,
* matrices, etc.
*
*/
struct BlockData {
// field
const string fieldName;
const string feName;
const string feNaturalBCName;
// material parameters
double mU; ///< magnetic constant N / A2
double ePsilon; ///< regularization paramater
// Natural boundary conditions
// Essential boundary conditions
int oRder; ///< approximation order
// Petsc data
DM dM;
Mat A;
Vec D, F;
: fieldName("MAGNETIC_POTENTIAL"), feName("MAGNETIC"),
feNaturalBCName("MAGENTIC_NATURAL_BC"), mU(1), ePsilon(0.1) {}
};
/**
* \brief get natural boundary conditions
* \ingroup maxwell_element
* @return error code
*/
if (bit->getName().compare(0, 9, "NATURALBC") == 0) {
Range faces;
CHKERR mField.get_moab().get_entities_by_type(bit->meshset, MBTRI,
faces, true);
faces, 1, true, blockData.naturalBc, moab::Interface::UNION);
blockData.naturalBc.merge(faces);
}
}
}
/**
* \brief get essential boundary conditions (only homogenous case is
* considered) \ingroup maxwell_element
* @return error code
*/
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
if (bit->getName().compare(0, 10, "ESSENTIALBC") == 0) {
Range faces;
CHKERR mField.get_moab().get_entities_by_type(bit->meshset, MBTRI,
faces, true);
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);
proc_skin, 1, true, blockData.essentialBc, moab::Interface::UNION);
blockData.essentialBc.merge(proc_skin);
}
}
/**
* \brief build problem data structures
* \ingroup maxwell_element
* @return error code
*/
// 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));
1);
3);
// meshset consisting all entities in mesh
EntityHandle root_set = mField.get_moab().get_root_set();
// // 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
// Set geometry approximation ordered
"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
// 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);
}
/**
* \brief create finite elements
* \ingroup maxwell_element
*
* Create volume and surface element. Surface element is used to integrate
* natural boundary conditions.
*
* @return error code
*/
// //Elements
"MESH_NODE_POSITIONS");
blockData.feNaturalBCName, "MESH_NODE_POSITIONS");
// build finite elemnts
}
/**
* \brief create problem
* \ingroup maxwell_element
*
* 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.
*
* @return error code
*/
// set up DM
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 DMSetUp(blockData.dM);
// remove essential DOFs
const MoFEM::Problem *problem_ptr;
CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
}
/**
* \brief destroy Distributed mesh manager
* \ingroup maxwell_element
* @return [description]
*/
CHKERR DMDestroy(&blockData.dM);
}
/** \brief solve problem
* \ingroup maxwell_element
*
* Create matrices; integrate over elements; solve linear system of equations
*
*/
MoFEMErrorCode solveProblem(const bool regression_test = false) {
VolumeFE vol_fe(mField);
auto material_det_vec = boost::make_shared<VectorDouble>();
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 mField.getInterface<MatrixManager>()
->createMPIAIJ<PetscGlobalIdx_mi_tag>(problem_ptr->getName(),
CHKERR MatZeroEntries(blockData.A);
CHKERR VecZeroEntries(blockData.F);
&vol_fe);
blockData.feNaturalBCName.c_str(), &tri_fe);
CHKERR MatAssemblyBegin(blockData.A, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(blockData.A, MAT_FINAL_ASSEMBLY);
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);
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);
}
/**
* \brief post-process results, i.e. save solution on the mesh
* \ingroup maxwell_element
* @return [description]
*/
PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(
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()));
&post_proc);
CHKERR post_proc.writeFile("out_values.h5m");
}
/** \brief calculate and assemble CurlCurl matrix
* \ingroup maxwell_element
\f[
\mathbf{A} = \int_\Omega \mu^{-1} \left( \nabla \times \mathbf{u} \cdot
\nabla \times \mathbf{v} \right) \textrm{d}\Omega \f] where \f[ \mathbf{u} =
\nabla \times \mathbf{B} \f] where \f$\mathbf{B}\f$ is magnetic flux and
\f$\mu\f$ is magnetic permeability.
For more details pleas look to \cite ivanyshyn2013computation
*/
struct OpCurlCurl
BlockData &blockData;
OpCurlCurl(BlockData &data)
: MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
data.fieldName, UserDataOperator::OPROWCOL),
blockData(data) {
sYmm = true;
}
MatrixDouble entityLocalMatrix;
/**
* \brief integrate matrix
* @param row_side local number of entity on element for row of the matrix
* @param col_side local number of entity on element for col of the matrix
* @param row_type type of row entity (EDGE/TRIANGLE/TETRAHEDRON)
* @param col_type type of col entity (EDGE/TRIANGLE/TETRAHEDRON)
* @param row_data structure of data, like base functions and associated
* methods to access those data on rows
* @param col_data structure of data, like base functions and associated
* methods to access those data on rows
* @return error code
*/
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
if (row_type == MBVERTEX)
if (col_type == MBVERTEX)
const int nb_row_dofs = row_data.getN().size2() / 3;
if (nb_row_dofs == 0)
const int nb_col_dofs = col_data.getN().size2() / 3;
if (nb_col_dofs == 0)
entityLocalMatrix.resize(nb_row_dofs, nb_col_dofs, false);
entityLocalMatrix.clear();
if (nb_row_dofs != static_cast<int>(row_data.getFieldData().size())) {
SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Number of base functions and DOFs on entity is different on "
"rows %d!=%d",
nb_row_dofs, row_data.getFieldData().size());
}
if (nb_col_dofs != static_cast<int>(col_data.getFieldData().size())) {
SETERRQ2(
PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Number of base functions and DOFs on entity is different on cols",
nb_col_dofs, col_data.getFieldData().size());
}
MatrixDouble row_curl_mat, col_curl_mat;
const double c0 = 1. / blockData.mU;
const int nb_gauss_pts = row_data.getN().size1();
auto t_row_curl_base = row_data.getFTensor2DiffN<3, 3>();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
// get integration weight scaled by volume
double w = getGaussPts()(3, gg) * getVolume();
for (int aa = 0; aa != nb_row_dofs; aa++) {
t_row_curl(i) = levi_civita(j, i, k) * t_row_curl_base(j, k);
FTensor::Tensor0<double *> t_local_mat(&entityLocalMatrix(aa, 0), 1);
auto t_col_curl_base = col_data.getFTensor2DiffN<3, 3>(gg, 0);
for (int bb = 0; bb != nb_col_dofs; bb++) {
t_col_curl(i) = levi_civita(j, i, k) * t_col_curl_base(j, k);
t_local_mat += c0 * w * t_row_curl(i) * t_col_curl(i);
++t_local_mat;
++t_col_curl_base;
}
++t_row_curl_base;
}
}
CHKERR MatSetValues(blockData.A, nb_row_dofs, &row_data.getIndices()[0],
nb_col_dofs, &col_data.getIndices()[0],
if (row_side != col_side || row_type != col_type) {
entityLocalMatrix = trans(entityLocalMatrix);
CHKERR MatSetValues(blockData.A, nb_col_dofs, &col_data.getIndices()[0],
nb_row_dofs, &row_data.getIndices()[0],
}
}
};
/** \brief calculate and assemble stabilization matrix
* \ingroup maxwell_element
\f[
\mathbf{A} = \int_\Omega \epsilon \mathbf{u} \cdot \mathbf{v}
\textrm{d}\Omega \f] where \f$\epsilon\f$ is regularization parameter.
*/
struct OpStab
BlockData &blockData;
OpStab(BlockData &data)
: MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
data.fieldName, UserDataOperator::OPROWCOL),
blockData(data) {
sYmm = true;
}
MatrixDouble entityLocalMatrix;
/**
* \brief integrate matrix
* @param row_side local number of entity on element for row of the matrix
* @param col_side local number of entity on element for col of the matrix
* @param row_type type of row entity (EDGE/TRIANGLE/TETRAHEDRON)
* @param col_type type of col entity (EDGE/TRIANGLE/TETRAHEDRON)
* @param row_data structure of data, like base functions and associated
* methods to access those data on rows
* @param col_data structure of data, like base functions and associated
* methods to access those data on rows
* @return error code
*/
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
if (row_type == MBVERTEX)
if (col_type == MBVERTEX)
const int nb_row_dofs = row_data.getN().size2() / 3;
if (nb_row_dofs == 0)
const int nb_col_dofs = col_data.getN().size2() / 3;
if (nb_col_dofs == 0)
entityLocalMatrix.resize(nb_row_dofs, nb_col_dofs, false);
entityLocalMatrix.clear();
if (nb_row_dofs != static_cast<int>(row_data.getFieldData().size())) {
SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Number of base functions and DOFs on entity is different on "
"rows %d!=%d",
nb_row_dofs, row_data.getFieldData().size());
}
if (nb_col_dofs != static_cast<int>(col_data.getFieldData().size())) {
SETERRQ2(
PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Number of base functions and DOFs on entity is different on cols",
nb_col_dofs, col_data.getFieldData().size());
}
MatrixDouble row_curl_mat, col_curl_mat;
const double c0 = 1. / blockData.mU;
const double c1 = blockData.ePsilon * c0;
const int nb_gauss_pts = row_data.getN().size1();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
// get integration weight scaled by volume
double w = getGaussPts()(3, gg) * getVolume();
&row_data.getVectorN<3>(gg)(0, HVEC0),
&row_data.getVectorN<3>(gg)(0, HVEC1),
&row_data.getVectorN<3>(gg)(0, HVEC2), 3);
for (int aa = 0; aa != nb_row_dofs; aa++) {
FTensor::Tensor0<double *> t_local_mat(&entityLocalMatrix(aa, 0), 1);
&col_data.getVectorN<3>(gg)(0, HVEC0),
&col_data.getVectorN<3>(gg)(0, HVEC1),
&col_data.getVectorN<3>(gg)(0, HVEC2), 3);
for (int bb = 0; bb != nb_col_dofs; bb++) {
t_local_mat += c1 * w * t_row_base(i) * t_col_base(i);
++t_col_base;
++t_local_mat;
}
++t_row_base;
}
}
// cerr << entityLocalMatrix << endl;
// cerr << endl;
CHKERR MatSetValues(blockData.A, nb_row_dofs, &row_data.getIndices()[0],
nb_col_dofs, &col_data.getIndices()[0],
if (row_side != col_side || row_type != col_type) {
entityLocalMatrix = trans(entityLocalMatrix);
CHKERR MatSetValues(blockData.A, nb_col_dofs, &col_data.getIndices()[0],
nb_row_dofs, &row_data.getIndices()[0],
}
}
};
/** \brief calculate essential boundary conditions
* \ingroup maxwell_element
\f[
\mathbf{F} = \int_{\partial\Omega} \ \mathbf{u} \cdot \mathbf{j}_i
\textrm{d}{\partial\Omega} \f] where \f$\mathbf{j}_i\f$ is current density
function.
Here simple current on coil is hard coded. In future more general
implementation is needed.
*/
struct OpNaturalBC
BlockData &blockData;
OpNaturalBC(BlockData &data)
data.fieldName, UserDataOperator::OPROW),
blockData(data) {}
VectorDouble naturalBC;
/**
* \brief integrate matrix
* \ingroup maxwell_element
* @param row_side local number of entity on element for row of the matrix
* @param row_type type of row entity (EDGE/TRIANGLE/TETRAHEDRON)
* @param row_data structure of data, like base functions and associated
* methods to access those data on rows
* @return error code
*/
MoFEMErrorCode doWork(int row_side, EntityType row_type,
if (row_type == MBVERTEX)
const int nb_row_dofs = row_data.getN().size2() / 3;
if (nb_row_dofs == 0)
naturalBC.resize(nb_row_dofs, false);
naturalBC.clear();
const int nb_gauss_pts = row_data.getN().size1();
auto t_row_base = row_data.getFTensor1N<3>();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
// get integration weight scaled by volume
double area;
area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
double w = getGaussPts()(2, gg) * area;
// Current is on surface where natural bc are applied. It is set that
// current is in XY plane, circular, around the coil.
const double x = getCoordsAtGaussPts()(gg, 0);
const double y = getCoordsAtGaussPts()(gg, 1);
const double r = sqrt(x * x + y * y);
t_j(0) = -y / r;
t_j(1) = +x / r;
t_j(2) = 0;
// double a = t_j(i) * t_tangent1(i);
// double b = t_j(i) * t_tangent2(i);
// t_j(i) = a * t_tangent1(i) + b * t_tangent2(i);
// ++t_tangent1;
// ++t_tangent2;
FTensor::Tensor0<double *> t_f(&naturalBC[0]);
for (int aa = 0; aa != nb_row_dofs; aa++) {
t_f += w * t_row_base(i) * t_j(i);
++t_row_base;
++t_f;
}
}
CHKERR VecSetValues(blockData.F, row_data.getIndices().size(),
}
};
/** \brief calculate and assemble CurlCurl matrix
* \ingroup maxwell_element
*/
struct OpPostProcessCurl
BlockData &blockData;
moab::Interface &postProcMesh;
std::vector<EntityHandle> &mapGaussPts;
OpPostProcessCurl(BlockData &data, moab::Interface &post_proc_mesh,
std::vector<EntityHandle> &map_gauss_pts)
: MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
data.fieldName, UserDataOperator::OPROW),
blockData(data), postProcMesh(post_proc_mesh),
mapGaussPts(map_gauss_pts) {}
MoFEMErrorCode doWork(int row_side, EntityType row_type,
if (row_type == MBVERTEX)
Tag th;
double def_val[] = {0, 0, 0};
CHKERR postProcMesh.tag_get_handle("MAGNETIC_INDUCTION_FIELD", 3,
MB_TYPE_DOUBLE, th,
MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
const int nb_row_dofs = row_data.getN().size2() / 3;
if (nb_row_dofs == 0)
const void *tags_ptr[mapGaussPts.size()];
if (nb_row_dofs != row_data.getFieldData().size())
SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong number of base functions and DOFs %d != %d",
nb_row_dofs, row_data.getFieldData().size());
const int nb_gauss_pts = row_data.getN().size1();
if (nb_gauss_pts != static_cast<int>(mapGaussPts.size())) {
SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Inconsistency number of dofs %d!=%d", nb_gauss_pts,
mapGaussPts.size());
}
CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0],
mapGaussPts.size(), tags_ptr);
auto t_curl_base = row_data.getFTensor2DiffN<3, 3>();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
// get pointer to tag values on entity (i.e. vertex on refined
// post-processing mesh)
double *ptr = &((double *)tags_ptr[gg])[0];
&ptr[2]);
// calculate curl value
auto t_dof = row_data.getFTensor0FieldData();
for (int aa = 0; aa != nb_row_dofs; aa++) {
t_curl(i) += t_dof * (levi_civita(j, i, k) * t_curl_base(j, k));
++t_curl_base;
++t_dof;
}
}
}
};
};
#endif //__MAGNETICELEMENT_HPP__
/******************************************************************************
* \defgroup maxwell_element Magnetic/Maxwell element
* \ingroup user_modules
******************************************************************************/
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MagneticElement::BlockData::~BlockData
~BlockData()
Definition: MagneticElement.hpp:91
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1631
FTensor::Tensor1< double, 3 >
EntityHandle
MagneticElement::BlockData::D
Vec D
Definition: MagneticElement.hpp:86
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Definition: HODataOperators.hpp:764
MagneticElement::BlockData::mU
double mU
magnetic constant N / A2
Definition: MagneticElement.hpp:72
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
MagneticElement::BlockData::dM
DM dM
Definition: MagneticElement.hpp:84
help
static char help[]
Definition: activate_deactivate_dofs.cpp:13
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MagneticElement::destroyProblem
MoFEMErrorCode destroyProblem()
destroy Distributed mesh manager
Definition: MagneticElement.hpp:307
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
MagneticElement::BlockData::oRder
int oRder
approximation order
Definition: MagneticElement.hpp:81
TriFE
Definition: prism_elements_from_surface.cpp:91
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
BasicFiniteElements.hpp
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator::getVolume
double getVolume() const
element volume (linear geometry)
Definition: VolumeElementForcesAndSourcesCore.hpp:161
MagneticElement::BlockData::essentialBc
Range essentialBc
Definition: MagneticElement.hpp:79
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition: Levi_Civita.hpp:617
MoFEM::DataOperator::doWork
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: DataOperators.hpp:40
MagneticElement::MagneticElement
MagneticElement(MoFEM::Interface &m_field)
Definition: MagneticElement.hpp:56
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1576
HVEC1
@ HVEC1
Definition: definitions.h:186
sdf.r
int r
Definition: sdf.py:8
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::ForcesAndSourcesCore::UserDataOperator::getGaussPts
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition: ForcesAndSourcesCore.hpp:1235
MagneticElement::BlockData::fieldName
const string fieldName
Definition: MagneticElement.hpp:67
MoFEM::VolumeElementForcesAndSourcesCore::VolumeElementForcesAndSourcesCore
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)
Definition: VolumeElementForcesAndSourcesCore.cpp:9
MagneticElement::BlockData::BlockData
BlockData()
Definition: MagneticElement.hpp:88
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
MagneticElement::BlockData::F
Vec F
Definition: MagneticElement.hpp:86
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator::getNormalsAtGaussPts
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
Definition: FaceElementForcesAndSourcesCore.hpp:290
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MagneticElement::getEssentialBc
MoFEMErrorCode getEssentialBc()
get essential boundary conditions (only homogenous case is considered)
Definition: MagneticElement.hpp:121
MagneticElement::createFields
MoFEMErrorCode createFields()
build problem data structures
Definition: MagneticElement.hpp:158
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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
BlockData
Definition: gel_analysis.cpp:74
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
MagneticElement::VolumeFE::VolumeFE
VolumeFE(MoFEM::Interface &m_field)
Definition: MagneticElement.hpp:34
MagneticElement::BlockData::naturalBc
Range naturalBc
Definition: MagneticElement.hpp:76
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
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
MagneticElement::TriFE::TriFE
TriFE(MoFEM::Interface &m_field)
Definition: MagneticElement.hpp:51
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MagneticElement::BlockData::feNaturalBCName
const string feNaturalBCName
Definition: MagneticElement.hpp:69
MoFEM::ForcesAndSourcesCore::UserDataOperator::getCoordsAtGaussPts
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
Definition: ForcesAndSourcesCore.hpp:1264
MagneticElement::postProcessResults
MoFEMErrorCode postProcessResults()
post-process results, i.e. save solution on the mesh
Definition: MagneticElement.hpp:406
MoFEM::DMMoFEMCreateMoFEM
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:114
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
main
int main(int argc, char *argv[])
Definition: activate_deactivate_dofs.cpp:15
FaceElementForcesAndSourcesCore
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
FTensor::Index< 'i', 3 >
MagneticElement::TriFE::getRule
int getRule(int order)
Definition: MagneticElement.hpp:53
MagneticElement::blockData
BlockData blockData
Definition: MagneticElement.hpp:94
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
MagneticElement::createProblem
MoFEMErrorCode createProblem()
create problem
Definition: MagneticElement.hpp:278
Range
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MagneticElement::mField
MoFEM::Interface & mField
Definition: MagneticElement.hpp:30
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
FTensor::Tensor0
Definition: Tensor0.hpp:16
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:71
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
MagneticElement.hpp
Implementation of magnetic element.
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MagneticElement::BlockData::ePsilon
double ePsilon
regularization paramater
Definition: MagneticElement.hpp:73
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MagneticElement::getNaturalBc
MoFEMErrorCode getNaturalBc()
get natural boundary conditions
Definition: MagneticElement.hpp:101
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:426
MagneticElement::BlockData::feName
const string feName
Definition: MagneticElement.hpp:68
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
sdf_wavy_2d.w
int w
Definition: sdf_wavy_2d.py:7
MoFEM::Types::BitRefLevel
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
MagneticElement::BlockData::A
Mat A
Definition: MagneticElement.hpp:85
MagneticElement::createElements
MoFEMErrorCode createElements()
create finite elements
Definition: MagneticElement.hpp:234
MagneticElement
Implementation of magneto-static problem (basic Implementation)
Definition: MagneticElement.hpp:28
MoFEM::Problem
Definition: ProblemsMultiIndices.hpp:54
HVEC2
@ HVEC2
Definition: definitions.h:186
MoFEM::CoreInterface::set_field_order
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
MoFEM::Problem::getName
auto getName() const
Definition: ProblemsMultiIndices.hpp:372
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MagneticElement::VolumeFE::getRule
int getRule(int order)
Definition: MagneticElement.hpp:36
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
tol
double tol
Definition: mesh_smoothing.cpp:27
HVEC0
@ HVEC0
Definition: definitions.h:186
MagneticElement::solveProblem
MoFEMErrorCode solveProblem(const bool regression_test=false)
solve problem
Definition: MagneticElement.hpp:319
MagneticElement::~MagneticElement
virtual ~MagneticElement()=default
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698