v0.13.1
Loading...
Searching...
No Matches
MagneticElement.hpp
/** \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 {
VolumeFE(MoFEM::Interface &m_field)
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
*/
MoFEMErrorCode getNaturalBc() {
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
*/
MoFEMErrorCode getEssentialBc() {
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
*/
MoFEMErrorCode createFields() {
// 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
*/
MoFEMErrorCode createElements() {
// //Elements
"MESH_NODE_POSITIONS");
blockData.feNaturalBCName, "MESH_NODE_POSITIONS");
// build finite elemnts
// build adjacencies
CHKERR mField.build_adjacencies(BitRefLevel().set(0));
}
/**
* \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
*/
MoFEMErrorCode createProblem() {
// 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.c_str());
CHKERR DMMoFEMAddElement(blockData.dM, blockData.feNaturalBCName.c_str());
CHKERR DMSetUp(blockData.dM);
// remove essential DOFs
const MoFEM::Problem *problem_ptr;
CHKERR DMMoFEMGetProblemPtr(blockData.dM, &problem_ptr);
CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
}
/**
* \brief destroy Distributed mesh manager
* \ingroup maxwell_element
* @return [description]
*/
MoFEMErrorCode destroyProblem() {
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 DMMoFEMGetProblemPtr(blockData.dM, &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);
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);
}
/**
* \brief post-process results, i.e. save solution on the mesh
* \ingroup maxwell_element
* @return [description]
*/
MoFEMErrorCode postProcessResults() {
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()));
CHKERR DMoFEMLoopFiniteElements(blockData.dM, blockData.feName.c_str(),
&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,
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
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],
&entityLocalMatrix(0, 0), ADD_VALUES);
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],
&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
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,
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
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],
&entityLocalMatrix(0, 0), ADD_VALUES);
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],
&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
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,
EntitiesFieldData::EntData &row_data) {
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(),
&row_data.getIndices()[0], &naturalBC[0], ADD_VALUES);
}
};
/** \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,
EntitiesFieldData::EntData &row_data) {
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
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BLOCKSET
Definition: definitions.h:148
@ HVEC0
Definition: definitions.h:186
@ HVEC1
Definition: definitions.h:186
@ HVEC2
Definition: definitions.h:186
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
virtual MoFEMErrorCode 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.
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
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:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
double w(const double g, const double t)
const double r
rate factor
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
int oRder
approximation order
double mU
magnetic constant N / A2
double ePsilon
regularization paramater
Implementation of magneto-static problem (basic Implementation)
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
virtual ~MagneticElement()=default
MoFEMErrorCode solveProblem(const bool regression_test=false)
solve problem
virtual moab::Interface & get_moab()=0
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 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.
Deprecated interface functions.
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
keeps basic data about problem
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.