#ifndef __MAGNETICELEMENT_HPP__
#define __MAGNETICELEMENT_HPP__
};
};
};
if (
bit->getName().compare(0, 9,
"NATURALBC") == 0) {
faces, true);
}
}
}
ParallelComm *pcomm =
if (
bit->getName().compare(0, 10,
"ESSENTIALBC") == 0) {
faces, true);
}
}
CHKERR skin.find_skin(0, tets,
false, skin_faces);
CHKERR pcomm->filter_pstatus(skin_faces,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &proc_skin);
}
}
0, 3, BitRefLevel().set(0));
1);
3);
"MESH_NODE_POSITIONS");
Projection10NodeCoordsOnField ent_method_material(
mField,
"MESH_NODE_POSITIONS");
}
"MESH_NODE_POSITIONS");
}
CHKERR DMRegister_MoFEM(
"DMMOFEM");
BitRefLevel().set(0));
}
}
MoFEMErrorCode
solveProblem(
const bool regression_test =
false) {
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));
new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
new OpHOSetCovariantPiolaTransformOnFace3D(
HCURL));
->createMPIAIJ<PetscGlobalIdx_mi_tag>(problem_ptr->
getName(),
&vol_fe);
KSP solver;
CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
CHKERR KSPSetFromOptions(solver);
SCATTER_REVERSE);
if (regression_test) {
double nrm2;
const double expected_value = 4.6772e+01;
if ((nrm2 - expected_value) / expected_value >
tol) {
"Regression test failed %6.4e != %6.4e", nrm2, expected_value);
}
}
}
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(
using OpPPMap = OpPostProcMapInMoab<3, 3>;
post_proc.getOpPtrVector().push_back(
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");
}
struct OpCurlCurl
blockData(data) {
sYmm = true;
}
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())) {
"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(
"Number of base functions and DOFs on entity is different on cols",
nb_col_dofs, col_data.getFieldData().size());
}
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++) {
for (int aa = 0; aa != nb_row_dofs; aa++) {
auto t_col_curl_base = col_data.getFTensor2DiffN<3, 3>(gg, 0);
for (int bb = 0; bb != nb_col_dofs; bb++) {
t_local_mat += c0 *
w * t_row_curl(
i) * t_col_curl(
i);
++t_local_mat;
++t_col_curl_base;
}
++t_row_curl_base;
}
}
nb_col_dofs, &col_data.getIndices()[0],
&entityLocalMatrix(0, 0), ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
entityLocalMatrix = trans(entityLocalMatrix);
nb_row_dofs, &row_data.getIndices()[0],
&entityLocalMatrix(0, 0), ADD_VALUES);
}
}
};
struct OpStab
blockData(data) {
sYmm = true;
}
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())) {
"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(
"Number of base functions and DOFs on entity is different on cols",
nb_col_dofs, col_data.getFieldData().size());
}
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++) {
&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;
}
}
nb_col_dofs, &col_data.getIndices()[0],
&entityLocalMatrix(0, 0), ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
entityLocalMatrix = trans(entityLocalMatrix);
nb_row_dofs, &row_data.getIndices()[0],
&entityLocalMatrix(0, 0), ADD_VALUES);
}
}
};
struct OpNaturalBC
blockData(data) {}
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++) {
double area;
const double r = sqrt(x * x + y * y);
t_j(2) = 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;
}
}
&row_data.getIndices()[0], &naturalBC[0], ADD_VALUES);
}
};
struct OpPostProcessCurl
moab::Interface &postProcMesh;
std::vector<EntityHandle> &mapGaussPts;
OpPostProcessCurl(
BlockData &data, moab::Interface &post_proc_mesh,
std::vector<EntityHandle> &map_gauss_pts)
blockData(data), postProcMesh(post_proc_mesh),
mapGaussPts(map_gauss_pts) {}
EntitiesFieldData::EntData &row_data) {
if (row_type == MBVERTEX)
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())
"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())) {
"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++) {
double *ptr = &((double *)tags_ptr[gg])[0];
&ptr[2]);
auto t_dof = row_data.getFTensor0FieldData();
for (int aa = 0; aa != nb_row_dofs; aa++) {
++t_curl_base;
++t_dof;
}
}
}
};
};
#endif
ForcesAndSourcesCore::UserDataOperator UserDataOperator
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ HCURL
field with continuous tangents
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
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 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 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_row(const std::string &fe_name, const std::string name_row)=0
set field row 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.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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.
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
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.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
int oRder
approximation order
const string feNaturalBCName
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.
default operator for TRI element
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
std::string meshPositionsFieldName
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
boost::ptr_deque< 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.
double getVolume() const
element volume (linear geometry)
Volume finite element base.