v0.13.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
*
*/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#include <BasicFiniteElements.hpp>
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();
// Read mesh to MOAB
const char *option;
option = "PARALLEL=READ_PART;"
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION;";
CHKERR moab.load_file(mesh_file_name, 0, option);
// 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;
}
const std::string default_options
std::string param_file
Implementation of magnetic element.
int main(int argc, char *argv[])
static char help[]
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
#define CHKERR
Inline error check.
Definition: definitions.h:548
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
CoreTmp< 0 > Core
Definition: Core.hpp:1096
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
Implementation of magneto-static problem (basic Implementation)
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
/** \file MagneticElement.hpp
* \brief Implementation of magnetic element
* \ingroup maxwell_element
* \example MagneticElement.hpp
*
*/
/*
* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#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
Range naturalBc;
// Essential boundary conditions
Range essentialBc;
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);
CHKERR mField.get_moab().get_adjacencies(
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);
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);
}
}
/**
* \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));
// add fields
1);
3);
// meshset consisting all entities in mesh
EntityHandle root_set = mField.get_moab().get_root_set();
// add entities to field
// // 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
// build adjacencies
}
/**
* \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);
// add elements to 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_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 mField.getInterface<MatrixManager>()
->createMPIAIJ<PetscGlobalIdx_mi_tag>(problem_ptr->getName(),
CHKERR MatZeroEntries(blockData.A);
CHKERR VecZeroEntries(blockData.F);
CHKERR VecGhostUpdateBegin(blockData.F, ADD_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(blockData.F, ADD_VALUES, SCATTER_FORWARD);
&vol_fe);
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);
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]
*/
CHKERR post_proc.generateReferenceElementMesh();
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", post_proc, false, true, false, true);
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));
&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
data.fieldName, UserDataOperator::OPROWCOL),
blockData(data) {
sYmm = true;
}
/**
* \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);
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);
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],
&entityLocalMatrix(0, 0), ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
CHKERR MatSetValues(blockData.A, nb_col_dofs, &col_data.getIndices()[0],
nb_row_dofs, &row_data.getIndices()[0],
&entityLocalMatrix(0, 0), ADD_VALUES);
}
}
};
/** \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
data.fieldName, UserDataOperator::OPROWCOL),
blockData(data) {
sYmm = true;
}
/**
* \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);
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++) {
&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],
&entityLocalMatrix(0, 0), ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
CHKERR MatSetValues(blockData.A, nb_col_dofs, &col_data.getIndices()[0],
nb_row_dofs, &row_data.getIndices()[0],
&entityLocalMatrix(0, 0), ADD_VALUES);
}
}
};
/** \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
data.fieldName, UserDataOperator::OPROW),
blockData(data) {}
/**
* \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;
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(),
&row_data.getIndices()[0], &naturalBC[0], ADD_VALUES);
}
};
/** \brief calculate and assemble CurlCurl matrix
* \ingroup maxwell_element
*/
struct OpPostProcessCurl
std::vector<EntityHandle> &mapGaussPts;
std::vector<EntityHandle> &map_gauss_pts)
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
******************************************************************************/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
EntitiesFieldData::EntData EntData
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ H1
continuous field
Definition: definitions.h:98
@ HCURL
field with continuous tangents
Definition: definitions.h:99
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ BLOCKSET
Definition: definitions.h:161
@ HVEC0
Definition: definitions.h:199
@ HVEC1
Definition: definitions.h:199
@ HVEC2
Definition: definitions.h:199
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1061
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: DMMMoFEM.cpp:130
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:461
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:481
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:384
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:544
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
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 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_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
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.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
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.
#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
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
double tol
const FTensor::Tensor2< T, Dim, Dim > Vec
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
double w(const double g, const double t)
const double r
rate factor
int oRder
approximation order
double mU
magnetic constant N / A2
double ePsilon
regularization paramater
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
integrate matrix
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
integrate matrix
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
BlockData & blockData
OpPostProcessCurl(BlockData &data, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
integrate matrix
TriFE(MoFEM::Interface &m_field)
VolumeFE(MoFEM::Interface &m_field)
MoFEMErrorCode createFields()
build problem data structures
MoFEMErrorCode createElements()
create finite elements
MoFEMErrorCode destroyProblem()
destroy Distributed mesh manager
MoFEMErrorCode postProcessResults()
post-process results, i.e. save solution on the mesh
MoFEMErrorCode createProblem()
create problem
MoFEM::Interface & mField
MoFEMErrorCode getEssentialBc()
get essential boundary conditions (only homogenous case is considered)
MoFEMErrorCode getNaturalBc()
get natural boundary conditions
MagneticElement(MoFEM::Interface &m_field)
virtual ~MagneticElement()=default
MoFEMErrorCode solveProblem(const bool regression_test=false)
solve problem
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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.
virtual moab::Interface & get_moab()=0
bool sYmm
If true assume that matrix is symmetric structure.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
keeps basic data about problem
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Post processing.