#include <BasicFiniteElements.hpp>
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[]) {
"-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";
std::ofstream file(
param_file.c_str(), std::ios::ate);
if (file.is_open()) {
file.close();
}
}
try {
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());
PetscBool flg_file;
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",
CHKERR PetscOptionsInt(
"-my_order",
"default approximation order",
"",
"-regression_test",
"if set norm of solution vector is check agains expected value ",
"",
PETSC_FALSE, ®ression_test, PETSC_NULL);
ierr = PetscOptionsEnd();
const char *option;
option = "PARALLEL=READ_PART;"
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION;";
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();
CHKERR moab.write_file(
"solution.h5m",
"MOAB",
"PARALLEL=WRITE_PART");
}
return 0;
}
const std::string default_options
Implementation of magnetic element.
int main(int argc, char *argv[])
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
#define CHKERR
Inline error check.
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
implementation of Data Operators for Forces and Sources
DeprecatedCoreInterface Interface
Implementation of magneto-static problem (basic Implementation)
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
#ifndef __MAGNETICELEMENT_HPP__
#define __MAGNETICELEMENT_HPP__
};
};
};
if (
bit->getName().compare(0, 9,
"NATURALBC") == 0) {
Range faces;
faces, true);
}
}
}
ParallelComm *pcomm =
if (
bit->getName().compare(0, 10,
"ESSENTIALBC") == 0) {
Range faces;
faces, true);
}
}
Range tets;
Range skin_faces;
CHKERR skin.find_skin(0, tets,
false, skin_faces);
Range proc_skin;
CHKERR pcomm->filter_pstatus(skin_faces,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &proc_skin);
}
}
1);
3);
"MESH_NODE_POSITIONS");
Projection10NodeCoordsOnField ent_method_material(
mField,
"MESH_NODE_POSITIONS");
}
"MESH_NODE_POSITIONS");
}
}
}
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>();
vol_fe.getOpPtrVector().push_back(
new OpCurlCurl(
blockData));
vol_fe.getOpPtrVector().push_back(
new OpStab(
blockData));
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));
->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);
}
}
}
CHKERR post_proc.generateReferenceElementMesh();
CHKERR post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
post_proc.getOpPtrVector().push_back(new OpPostProcessCurl(
blockData, post_proc.postProcMesh, post_proc.mapGaussPts));
&post_proc);
CHKERR post_proc.writeFile(
"out_values.h5m");
}
struct OpCurlCurl
}
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)
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 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],
if (row_side != col_side || row_type != col_type) {
nb_row_dofs, &row_data.getIndices()[0],
}
}
};
struct OpStab
}
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)
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 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],
if (row_side != col_side || row_type != col_type) {
nb_row_dofs, &row_data.getIndices()[0],
}
}
};
struct OpNaturalBC
if (row_type == MBVERTEX)
const int nb_row_dofs = row_data.getN().size2() / 3;
if (nb_row_dofs == 0)
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;
area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
double w = getGaussPts()(2, gg) * area;
const double x = getCoordsAtGaussPts()(gg, 0);
const double y = getCoordsAtGaussPts()(gg, 1);
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
std::vector<EntityHandle> &map_gauss_pts)
if (row_type == MBVERTEX)
double def_val[] = {0, 0, 0};
MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
const int nb_row_dofs = row_data.getN().size2() / 3;
if (nb_row_dofs == 0)
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,
}
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
EntitiesFieldData::EntData EntData
@ 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 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()
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
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.
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
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.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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.
UBlasMatrix< double > MatrixDouble
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
UBlasVector< double > VectorDouble
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
constexpr auto VecSetValues
constexpr auto MatSetValues
const double r
rate factor
int oRder
approximation order
const string feNaturalBCName
double mU
magnetic constant N / A2
double ePsilon
regularization paramater
MatrixDouble entityLocalMatrix
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
integrate matrix
OpCurlCurl(BlockData &data)
OpNaturalBC(BlockData &data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
integrate matrix
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
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)
MatrixDouble entityLocalMatrix
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.
double getVolume() const
element volume (linear geometry)
Volume finite element with switches.