static char help[] =
"...\n\n";
};
};
using PostProcEle::PostProcEle;
protected:
};
private:
};
}
PetscBool load_file = PETSC_FALSE;
PETSC_NULL);
if (load_file == PETSC_FALSE) {
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);
for (auto d : {1, 2})
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);
CHKERR moab.get_adjacencies(&tri, 1, 1,
true, adj);
}
} else {
}
}
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)
}
}
auto post_proc_fe = boost::make_shared<MyPostProc>(
mField);
post_proc_fe->generateReferenceElementMesh();
auto jac_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
}
}
{
auto u_ptr = boost::make_shared<VectorDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"U", u_ptr}},
{},
{},
{}
)
);
} break;
{
auto u_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{},
{{"U", u_ptr}},
{},
{}
)
);
} break;
default:
break;
}
auto scale_tag_val = [&]() {
auto &post_proc_mesh = post_proc_fe->getPostProcMesh();
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 post_proc_fe->writeFile(
"out_base_dof_" + boost::lexical_cast<std::string>(nb) + ".h5m");
CHKERR post_proc_fe->getPostProcMesh().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(
LogManager::createSink(LogManager::getStrmWorld(), "PLOTBASE"));
LogManager::setLog("PLOTBASE");
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
}
return 0;
}
moab::Core core_ref;
moab::Interface &moab_ref = core_ref;
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,
"");
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);
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);
for (auto d : {0, 1, 2})
gaussPts(d, gg) = coords[
d];
nodes_pts_map[*nit] = gg;
}
}
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) {
}
}
}
const int num_nodes = gaussPts.size2();
switch (numeredEntFiniteElementPtr->getEntType()) {
case MBTRI:
&gaussPts(0, 0), &gaussPts(1, 0), num_nodes);
break;
case MBQUAD: {
for (int gg = 0; gg != num_nodes; gg++) {
double ksi = gaussPts(0, gg);
double eta = gaussPts(1, gg);
}
} break;
case MBTET: {
&gaussPts(0, 0), &gaussPts(1, 0),
&gaussPts(2, 0), num_nodes);
} break;
case MBHEX: {
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 getPostProcMesh().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 getPostProcMesh().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_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 getPostProcMesh().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]);
}
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;
}
}
ParallelComm *pcomm_post_proc_mesh =
if (pcomm_post_proc_mesh != NULL)
delete pcomm_post_proc_mesh;
};
auto resolve_shared_ents = [&]() {
ParallelComm *pcomm_post_proc_mesh =
if (pcomm_post_proc_mesh == NULL) {
pcomm_post_proc_mesh = new ParallelComm(
&(getPostProcMesh()),
PETSC_COMM_WORLD );
}
CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
};
}
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
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#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 ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'i', SPACE_DIM > i
const 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)
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)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
double zeta
Viscous hardening.
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
MoFEMErrorCode readMesh()
[Run problem]
MoFEMErrorCode setIntegrationRules()
[Set up problem]
FieldApproximationBase base
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode solveSystem()
[Solve]
MoFEMErrorCode createCommonData()
[Set up problem]
MoFEMErrorCode runProblem()
[Run problem]
MoFEM::Interface & mField
MoFEMErrorCode setupProblem()
[Run problem]
MoFEMErrorCode outputResults()
[Solve]
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode rebuild_database(int verb=DEFAULT_VERBOSITY)=0
Clear database and initialize it once again.
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)
Get vector field for H-div approximation.
Get value at integration points for scalar field.
Make Hdiv space from Hcurl space in 2d.
Post post-proc data at points from hash maps.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Simple interface for fast problem set-up.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
MoFEMErrorCode postProcess()
MoFEMErrorCode generateReferenceElementMesh()
ublas::matrix< int > refEleMap
MoFEMErrorCode setGaussPts(int order)
MoFEMErrorCode preProcess()
MatrixDouble shapeFunctions