static char help[] =
"...\n\n";
#include <BasicFiniteElements.hpp>
using DomainEle = FaceElementForcesAndSourcesCoreBase;
};
using DomainEle = VolumeElementForcesAndSourcesCoreBase;
};
protected:
ublas::matrix<int> refEleMap;
};
private:
Simple *simpleInterface;
};
}
PetscBool load_file = PETSC_FALSE;
PETSC_NULL);
if (load_file == PETSC_FALSE) {
auto &moab = mField.get_moab();
double tet_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
for (int nn = 0; nn < 4; nn++) {
CHKERR moab.create_vertex(&tet_coords[3 * nn], nodes[nn]);
}
CHKERR moab.create_element(MBTET, nodes, 4, tet);
Range adj;
CHKERR moab.get_adjacencies(&tet, 1,
d,
true, adj);
}
double tri_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
for (int nn = 0; nn < 3; nn++) {
CHKERR moab.create_vertex(&tri_coords[3 * nn], nodes[nn]);
}
CHKERR moab.create_element(MBTRI, nodes, 3, tri);
Range adj;
CHKERR moab.get_adjacencies(&tri, 1, 1,
true, adj);
}
CHKERR mField.rebuild_database();
CHKERR mField.getInterface(simpleInterface);
CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(
0,
SPACE_DIM, simpleInterface->getBitRefLevel());
} else {
CHKERR mField.getInterface(simpleInterface);
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
}
}
enum bases { AINSWORTH, AINSWORTH_LOBATTO, DEMKOWICZ, BERNSTEIN, LASBASETOP };
const char *list_bases[] = {"ainsworth", "ainsworth_lobatto", "demkowicz",
"bernstein"};
PetscBool flg;
PetscInt choice_base_value = AINSWORTH;
&choice_base_value, &flg);
if (flg != PETSC_TRUE)
if (choice_base_value == AINSWORTH)
if (choice_base_value == AINSWORTH_LOBATTO)
else if (choice_base_value == DEMKOWICZ)
else if (choice_base_value == BERNSTEIN)
enum spaces { H1SPACE, L2SPACE, HCURLSPACE, HDIVSPACE, LASBASETSPACE };
const char *list_spaces[] = {"h1", "l2", "hcurl", "hdiv"};
PetscInt choice_space_value = H1SPACE;
LASBASETSPACE, &choice_space_value, &flg);
if (flg != PETSC_TRUE)
if (choice_space_value == H1SPACE)
else if (choice_space_value == L2SPACE)
else if (choice_space_value == HCURLSPACE)
else if (choice_space_value == HDIVSPACE)
CHKERR simpleInterface->addDomainField(
"U", space, base, 1);
CHKERR simpleInterface->setUp();
}
}
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainLhsFE().reset();
auto post_proc_fe = boost::make_shared<MyPostProc>(mField);
post_proc_fe->generateReferenceElementMesh();
pipeline_mng->getDomainRhsFE() = post_proc_fe;
auto jac_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(new OpMakeHdivFromHcurl());
post_proc_fe->getOpPtrVector().push_back(
}
}
post_proc_fe->addFieldValuesPostProc("U");
auto scale_tag_val = [&]() {
auto &post_proc_mesh = post_proc_fe->postProcMesh;
Range nodes;
CHKERR post_proc_mesh.get_entities_by_type(0, MBVERTEX, nodes);
CHKERR post_proc_mesh.tag_get_handle(
"U",
th);
int length;
CHKERR post_proc_mesh.tag_get_length(
th, length);
std::vector<double> data(nodes.size() * length);
CHKERR post_proc_mesh.tag_get_data(
th, nodes, &*data.begin());
double max_v = 0;
for (
int i = 0;
i != nodes.size(); ++
i) {
for (
int d = 0;
d != length; ++
d)
v += pow(data[length *
i +
d], 2);
max_v = std::max(max_v,
v);
}
CHKERR post_proc_mesh.tag_set_data(
th, nodes, &*data.begin());
};
size_t nb = 0;
auto dofs_ptr = mField.get_dofs();
for (auto dof_ptr : (*dofs_ptr)) {
MOFEM_LOG(
"PLOTBASE", Sev::verbose) << *dof_ptr;
auto &val = const_cast<double &>(dof_ptr->getFieldData());
val = 1;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile(
"out_base_dof_" + boost::lexical_cast<std::string>(nb) + ".h5m");
CHKERR post_proc_fe->postProcMesh.delete_mesh();
val = 0;
++nb;
};
}
int main(
int argc,
char *argv[]) {
try {
DMType dm_name = "DMMOFEM";
auto core_log = logging::core::get();
core_log->add_sink(
}
return 0;
}
char ref_mesh_file_name[255];
strcpy(ref_mesh_file_name, "ref_mesh2d.h5m");
strcpy(ref_mesh_file_name, "ref_mesh3d.h5m");
else
"Dimension not implemented");
255, PETSC_NULL);
CHKERR moab_ref.load_file(ref_mesh_file_name, 0,
"");
Range elems;
CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
CHKERR moab_ref.add_entities(meshset, elems);
CHKERR moab_ref.convert_entities(meshset,
true,
false,
false);
CHKERR moab_ref.delete_entities(&meshset, 1);
Range elem_nodes;
CHKERR moab_ref.get_connectivity(elems, elem_nodes,
false);
std::map<EntityHandle, int> nodes_pts_map;
gaussPts.resize(
SPACE_DIM + 1, elem_nodes.size(),
false);
gaussPts.clear();
Range::iterator nit = elem_nodes.begin();
for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
double coords[3];
CHKERR moab_ref.get_coords(&*nit, 1, coords);
gaussPts(
d, gg) = coords[
d];
nodes_pts_map[*nit] = gg;
}
refEleMap.resize(elems.size(), 3 + 3);
refEleMap.resize(elems.size(), 4 + 6);
}
Range::iterator tit = elems.begin();
for (int tt = 0; tit != elems.end(); ++tit, ++tt) {
int num_nodes;
CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
for (int nn = 0; nn != num_nodes; ++nn) {
refEleMap(tt, nn) = nodes_pts_map[conn[nn]];
}
}
}
const int num_nodes = gaussPts.size2();
switch (numeredEntFiniteElementPtr->getEntType()) {
case MBTRI:
shapeFunctions.resize(num_nodes, 3);
&gaussPts(0, 0), &gaussPts(1, 0), num_nodes);
break;
case MBQUAD: {
shapeFunctions.resize(num_nodes, 4);
for (int gg = 0; gg != num_nodes; gg++) {
double ksi = gaussPts(0, gg);
double eta = gaussPts(1, gg);
}
} break;
case MBTET: {
shapeFunctions.resize(num_nodes, 8);
&gaussPts(0, 0), &gaussPts(1, 0),
&gaussPts(2, 0), num_nodes);
} break;
case MBHEX: {
shapeFunctions.resize(num_nodes, 8);
for (int gg = 0; gg != num_nodes; gg++) {
double ksi = gaussPts(0, gg);
double eta = gaussPts(1, gg);
double zeta = gaussPts(2, gg);
}
} break;
default:
"Not implemented element type");
}
ReadUtilIface *iface;
CHKERR postProcMesh.query_interface(iface);
std::vector<double *> arrays;
CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays);
mapGaussPts.resize(gaussPts.size2());
for (int gg = 0; gg != num_nodes; ++gg)
mapGaussPts[gg] = startv + gg;
int def_in_the_loop = -1;
CHKERR postProcMesh.tag_get_handle(
"NB_IN_THE_LOOP", 1, MB_TYPE_INTEGER,
th,
MB_TAG_CREAT | MB_TAG_SPARSE,
&def_in_the_loop);
const int num_el = refEleMap.size1();
const int num_nodes_on_ele = refEleMap.size2();
CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTRI, 0,
starte, conn);
CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTET, 0,
starte, conn);
else
"Dimension not implemented");
for (unsigned int tt = 0; tt != refEleMap.size1(); ++tt) {
for (int nn = 0; nn != num_nodes_on_ele; ++nn)
conn[num_nodes_on_ele * tt + nn] = mapGaussPts[refEleMap(tt, nn)];
}
CHKERR iface->update_adjacencies(starte, num_el, num_nodes_on_ele, conn);
auto physical_elements = Range(starte, starte + num_el - 1);
CHKERR postProcMesh.tag_clear_data(
th, physical_elements, &(nInTheLoop));
int fe_num_nodes;
{
mField.get_moab().get_connectivity(fe_ent, conn, fe_num_nodes, true);
coords.resize(3 * fe_num_nodes, false);
CHKERR mField.get_moab().get_coords(conn, fe_num_nodes, &coords[0]);
}
&*shapeFunctions.data().begin());
arrays[0], arrays[1], arrays[2]);
const double *t_coords_ele_x = &coords[0];
const double *t_coords_ele_y = &coords[1];
const double *t_coords_ele_z = &coords[2];
for (int gg = 0; gg != num_nodes; ++gg) {
t_coords_ele_x, t_coords_ele_y, t_coords_ele_z);
for (int nn = 0; nn != fe_num_nodes; ++nn) {
t_coords(
i) += t_n * t_ele_coords(
i);
for (auto ii : {0, 1, 2})
if (std::abs(t_coords(ii)) < std::numeric_limits<float>::epsilon())
t_coords(ii) = 0;
++t_ele_coords;
++t_n;
}
++t_coords;
}
}
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
#define EXECUTABLE_DIMENSION
int main(int argc, char *argv[])
[Check results]
EntitiesFieldData::EntData EntData
ElementsAndOps< SPACE_DIM >::PostProcEle PostProcEle
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode assembleSystem()
[Boundary condition]
MoFEMErrorCode readMesh()
[Run problem]
MoFEMErrorCode setIntegrationRules()
[Set up problem]
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode solveSystem()
[Solve]
MoFEMErrorCode createCommonData()
[Create common data]
MoFEMErrorCode runProblem()
[Create common data]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode outputResults()
[Solve]
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.
Data on single entity (This is passed as argument to DataOperator::doWork)
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
MoFEMErrorCode generateReferenceElementMesh()
MoFEMErrorCode setGaussPts(int order)