v0.13.1
Searching...
No Matches
plot_base.cpp

Utility for plotting base functions for different spaces, polynomial bases

/**
* \file plot_base.cpp
* \example plot_base.cpp
*
* Utility for plotting base functions for different spaces, polynomial bases
*/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
};
template <> struct ElementsAndOps<3> {
};
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
struct MyPostProc : public PostProcEle {
using PostProcEle::PostProcEle;
protected:
ublas::matrix<int> refEleMap;
};
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
private:
};
//! [Run programme]
}
//! [Run programme]
PETSC_NULL);
auto &moab = mField.get_moab();
if (SPACE_DIM == 3) {
// create one tet
double tet_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
EntityHandle nodes[4];
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})
}
if (SPACE_DIM == 2) {
// create one triangle
double tri_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
EntityHandle nodes[3];
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);
}
// Add all elements to database
CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(
0, SPACE_DIM, simpleInterface->getBitRefLevel());
} else {
CHKERR simpleInterface->getOptions();
}
}
//! [Set up problem]
// Declare elements
enum bases { AINSWORTH, AINSWORTH_LOBATTO, DEMKOWICZ, BERNSTEIN, LASBASETOP };
const char *list_bases[] = {"ainsworth", "ainsworth_lobatto", "demkowicz",
"bernstein"};
PetscBool flg;
PetscInt choice_base_value = AINSWORTH;
CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases, LASBASETOP,
&choice_base_value, &flg);
if (flg != PETSC_TRUE)
SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "base not set");
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;
CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-space", list_spaces,
LASBASETSPACE, &choice_space_value, &flg);
if (flg != PETSC_TRUE)
SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "space not set");
space = H1;
if (choice_space_value == H1SPACE)
space = H1;
else if (choice_space_value == L2SPACE)
space = L2;
else if (choice_space_value == HCURLSPACE)
else if (choice_space_value == HDIVSPACE)
int order = 3;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR simpleInterface->setFieldOrder("U", order);
}
//! [Set up problem]
//! [Set integration rule]
}
//! [Set integration rule]
//! [Create common data]
//! [Create common data]
//! [Boundary condition]
//! [Boundary condition]
//! [Push operators to pipeline]
//! [Push operators to pipeline]
//! [Solve]
//! [Solve]
pipeline_mng->getDomainLhsFE().reset();
auto post_proc_fe = boost::make_shared<MyPostProc>(mField);
post_proc_fe->generateReferenceElementMesh();
pipeline_mng->getDomainRhsFE() = post_proc_fe;
if (SPACE_DIM == 2) {
if (space == HCURL) {
auto jac_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateHOJacForFace(jac_ptr));
post_proc_fe->getOpPtrVector().push_back(new OpMakeHdivFromHcurl());
post_proc_fe->getOpPtrVector().push_back(
}
}
switch (space) {
case H1:
case L2:
{
auto u_ptr = boost::make_shared<VectorDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("U", u_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"U", u_ptr}},
{},
{},
{}
)
);
} break;
case HCURL:
case HDIV:
{
auto u_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
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();
Range nodes;
CHKERR post_proc_mesh.get_entities_by_type(0, MBVERTEX, nodes);
Tag th;
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) {
double v = 0;
for (int d = 0; d != length; ++d)
v += pow(data[length * i + d], 2);
v = std::sqrt(v);
max_v = std::max(max_v, v);
}
for (auto &v : data)
v /= max_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 scale_tag_val();
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;
};
}
//! [Postprocess results]
//! [Check results]
//! [Check results]
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
// Add logging channel for example
auto core_log = logging::core::get();
LogManager::createSink(LogManager::getStrmWorld(), "PLOTBASE"));
LogManager::setLog("PLOTBASE");
MOFEM_LOG_TAG("PLOTBASE", "plotbase");
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database insterface
//! [Create MoFEM]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
}
return 0;
}
moab::Core core_ref;
moab::Interface &moab_ref = core_ref;
char ref_mesh_file_name[255];
if (SPACE_DIM == 2)
strcpy(ref_mesh_file_name, "ref_mesh2d.h5m");
else if (SPACE_DIM == 3)
strcpy(ref_mesh_file_name, "ref_mesh3d.h5m");
else
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"Dimension not implemented");
CHKERR PetscOptionsGetString(PETSC_NULL, "", "-ref_file", ref_mesh_file_name,
255, PETSC_NULL);
// Get elements
Range elems;
CHKERR moab_ref.get_entities_by_dimension(0, SPACE_DIM, elems);
EntityHandle meshset;
CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
CHKERR moab_ref.convert_entities(meshset, true, false, false);
CHKERR moab_ref.delete_entities(&meshset, 1);
// Get nodes on the mesh
Range elem_nodes;
CHKERR moab_ref.get_connectivity(elems, elem_nodes, false);
// Map node entity and Gauss pint number
std::map<EntityHandle, int> nodes_pts_map;
// Set gauss points coordinates from the reference mesh
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;
}
if (SPACE_DIM == 2) {
// Set size of adjacency matrix (note ho order nodes 3 nodes and 3 nodes on
// edges)
refEleMap.resize(elems.size(), 3 + 3);
} else if (SPACE_DIM == 3) {
refEleMap.resize(elems.size(), 4 + 6);
}
Range::iterator tit = elems.begin();
for (int tt = 0; tit != elems.end(); ++tit, ++tt) {
const EntityHandle *conn;
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();
// Calculate shape functions
switch (numeredEntFiniteElementPtr->getEntType()) {
case MBTRI:
shapeFunctions.resize(num_nodes, 3);
CHKERR Tools::shapeFunMBTRI(&*shapeFunctions.data().begin(),
&gaussPts(0, 0), &gaussPts(1, 0), num_nodes);
break;
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);
CHKERR Tools::shapeFunMBTET(&*shapeFunctions.data().begin(),
&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);
shapeFunctions(gg, 0) = N_MBHEX0(ksi, eta, zeta);
shapeFunctions(gg, 1) = N_MBHEX1(ksi, eta, zeta);
shapeFunctions(gg, 2) = N_MBHEX2(ksi, eta, zeta);
shapeFunctions(gg, 3) = N_MBHEX3(ksi, eta, zeta);
shapeFunctions(gg, 4) = N_MBHEX4(ksi, eta, zeta);
shapeFunctions(gg, 5) = N_MBHEX5(ksi, eta, zeta);
shapeFunctions(gg, 6) = N_MBHEX6(ksi, eta, zeta);
shapeFunctions(gg, 7) = N_MBHEX7(ksi, eta, zeta);
}
} break;
default:
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"Not implemented element type");
}
// Create physical nodes
// MoAB interface allowing for creating nodes and elements in the bulk
CHKERR postProcMesh.query_interface(iface);
std::vector<double *> arrays; /// pointers to memory allocated by MoAB for
/// storing X, Y, and Z coordinates
EntityHandle startv; // Starting handle for vertex
// Allocate memory for num_nodes, and return starting handle, and access to
// memort.
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;
Tag th;
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);
// Create physical elements
const int num_el = refEleMap.size1();
const int num_nodes_on_ele = refEleMap.size2();
EntityHandle starte; // Starting handle to first created element
// Create tris/tets in the bulk in MoAB database
if (SPACE_DIM == 2)
CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTRI, 0,
starte, conn);
else if (SPACE_DIM == 3)
CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTET, 0,
starte, conn);
else
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"Dimension not implemented");
// At this point elements (memory for elements) is allocated, at code bellow
// actual connectivity of elements is set.
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)];
}
// and elements are ready to use.
auto physical_elements = Range(starte, starte + num_el - 1);
CHKERR postProcMesh.tag_clear_data(th, physical_elements, &(nInTheLoop));
EntityHandle fe_ent = numeredEntFiniteElementPtr->getEnt();
int fe_num_nodes;
{
const EntityHandle *conn;
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]);
}
// Set physical coordinates to physical nodes
&*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);
t_coords(i) = 0;
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;
}
}
moab::Interface &moab = coreMesh;
ParallelComm *pcomm_post_proc_mesh =
ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
if (pcomm_post_proc_mesh != NULL)
delete pcomm_post_proc_mesh;
};
auto resolve_shared_ents = [&]() {
ParallelComm *pcomm_post_proc_mesh =
ParallelComm::get_pcomm(&(postProcMesh), MYPCOMM_INDEX);
if (pcomm_post_proc_mesh == NULL) {
// wrapRefMeshComm =
// boost::make_shared<WrapMPIComm>(T::mField.get_comm(), false);
pcomm_post_proc_mesh = new ParallelComm(
&(postProcMesh),
PETSC_COMM_WORLD /*(T::wrapRefMeshComm)->get_comm()*/);
}
CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
};
CHKERR resolve_shared_ents();
}
static char help[]
int main()
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
FieldSpace
approximation spaces
Definition: definitions.h:82
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#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
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:35
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#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
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
Definition: fem_tools.h:60
#define N_MBHEX7(x, y, z)
Definition: fem_tools.h:78
#define N_MBHEX3(x, y, z)
Definition: fem_tools.h:74
#define N_MBHEX5(x, y, z)
Definition: fem_tools.h:76
#define N_MBHEX4(x, y, z)
Definition: fem_tools.h:75
#define N_MBHEX0(x, y, z)
Definition: fem_tools.h:71
#define N_MBHEX6(x, y, z)
Definition: fem_tools.h:77
#define N_MBHEX2(x, y, z)
Definition: fem_tools.h:73
Definition: fem_tools.h:57
#define N_MBHEX1(x, y, z)
Definition: fem_tools.h:72
Definition: fem_tools.h:59
Definition: fem_tools.h:58
ElementsAndOps< SPACE_DIM >::PostProcEle PostProcEle
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
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)
Definition: dTensor0.hpp:27
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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
Definition: plastic.cpp:10
double zeta
Definition: plastic.cpp:114
double eta
Definition: plastic.cpp:115
constexpr int SPACE_DIM
Definition: plot_base.cpp:26
[Example]
Definition: plastic.cpp:139
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: helmholtz.cpp:106
MoFEMErrorCode assembleSystem()
[Applying essential BC]
Definition: helmholtz.cpp:154
[run problem]
Definition: helmholtz.cpp:73
MoFEMErrorCode setIntegrationRules()
[Set up problem]
Definition: plot_base.cpp:203
FieldApproximationBase base
Definition: plot_base.cpp:68
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: contact.cpp:659
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
Definition: helmholtz.cpp:235
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:259
Simple * simpleInterface
Definition: helmholtz.cpp:45
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:174
FieldSpace space
Definition: plot_base.cpp:69
MoFEM::Interface & mField
Definition: plastic.cpp:146
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:186
MoFEMErrorCode outputResults()
[Solve]
Definition: helmholtz.cpp:255
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode rebuild_database(int verb=DEFAULT_VERBOSITY)=0
Clear database and initialize it once again.
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
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.
Definition: Simple.hpp:26
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: plot_base.cpp:34
MoFEMErrorCode postProcess()
Definition: plot_base.cpp:622
MoFEMErrorCode generateReferenceElementMesh()
Definition: plot_base.cpp:397
ublas::matrix< int > refEleMap
Definition: plot_base.cpp:44
MoFEMErrorCode setGaussPts(int order)
Definition: plot_base.cpp:468
MoFEMErrorCode preProcess()
Definition: plot_base.cpp:612
MatrixDouble shapeFunctions
Definition: plot_base.cpp:45