Adding prims on the surface and checking conformity between quads, triangles and prism
static char help[] =
"...\n\n";
static constexpr
double delta =
static constexpr std::array<double, 3>
d3 = {0, 0, 0};
static constexpr std::array<double, 3>
d4 = {0, 0,
delta};
inline static double getArg(double x) {
};
int x, y, z;
: x(getArg(coords[0])), y(getArg(coords[1])), z(getArg(coords[2])),
};
typedef multi_index_container<
indexed_by<
hashed_unique<composite_key<
member<CoordsAndHandle, int, &CoordsAndHandle::y>,
member<CoordsAndHandle, int, &CoordsAndHandle::z>>>
>>
public:
std::vector<EntityHandle> nodeHandles;
};
int getRuleTrianglesOnly(
int order);
int getRuleThroughThickness(
int order);
private:
};
template <
typename OP>
struct Op :
public OP {
public:
std::vector<EntityHandle> nodeHandles;
};
int getRule(int order_row, int order_col, int order_data);
MoFEMErrorCode setGaussPts(
int order_row,
int order_col,
int order_data);
private:
};
int getRule(int order_row, int order_col, int order_data);
MoFEMErrorCode setGaussPts(
int order_row,
int order_col,
int order_data);
private:
std::array<Range, 3> &edgeBlocks;
};
int getRule(int order_row, int order_col, int order_data);
MoFEMErrorCode setGaussPts(
int order_row,
int order_col,
int order_data);
private:
std::array<Range, 3> &edgeBlocks;
};
int main(
int argc,
char *argv[]) {
try {
int rank;
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
PetscBool flg = PETSC_TRUE;
255, &flg);
if (flg != PETSC_TRUE)
"error -my_file (MESH FILE NEEDED)");
const char *option;
option = "";
CHKERR moab.get_entities_by_type(0, MBTRI, tris,
false);
CHKERR moab.get_connectivity(tris, tris_verts);
CHKERR moab.get_coords(tris_verts, &tri_coords(0, 0));
std::array<Range, 3> edge_block;
for (auto b : {1, 2, 3})
Range extrude_prisms = prisms;
extrude_prisms, false, add_prims_layer);
prisms_from_surface_interface->
setThickness(add_prims_layer,
d3.data(),
extrude_prisms = add_prims_layer;
prisms.merge(add_prims_layer);
add_prims_layer.clear();
}
CHKERR moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
CHKERR moab.add_entities(meshset, prisms);
CHKERR moab.get_connectivity(prisms, verts);
CHKERR moab.get_coords(verts, &coords(0, 0));
for (
size_t v = 0;
v != verts.size(); ++
v)
CHKERR moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER,
one_prism_meshset);
std::array<double, 18> one_prism_coords = {0, 0, 0, 1, 0, 0, 0, 1, 0,
0, 0, 1, 1, 0, 1, 0, 1, 1};
std::array<EntityHandle, 6> one_prism_nodes;
for (
int n = 0;
n != 6; ++
n)
CHKERR moab.create_vertex(&one_prism_coords[3 *
n], one_prism_nodes[
n]);
CHKERR m_field.
get_moab().create_element(MBPRISM, one_prism_nodes.data(), 6,
one_prism);
one_prism_range.insert(one_prism);
CHKERR moab.add_entities(one_prism_meshset, one_prism_range);
CHKERR moab.get_connectivity(one_prism_range, one_prism_verts);
CHKERR moab.add_entities(one_prism_meshset, one_prism_verts);
Range one_prism_adj_ents;
for (
int d = 1;
d != 3; ++
d)
CHKERR moab.get_adjacencies(one_prism_range,
d,
true, one_prism_adj_ents,
moab::Interface::UNION);
CHKERR moab.add_entities(one_prism_meshset, one_prism_adj_ents);
bit_level0.set(0);
one_prism_range, bit_level0);
bit_level0);
"EDGE");
"TRI");
"QUAD");
MBPRISM, "PRISM");
PrismFE fe_prism(m_field, tri_coords);
fe_prism.getOpPtrVector().push_back(
new PrismOp(moab, map_coords));
EdgeFE fe_edge(m_field, edge_block, one_prism);
fe_edge.getOpPtrVector().push_back(
moab, map_coords, one_prism));
TriFE fe_tri(m_field, tri_coords, one_prism);
fe_tri.getOpPtrVector().push_back(
moab, map_coords, one_prism));
QuadFE fe_quad(m_field, edge_block, one_prism);
fe_quad.getOpPtrVector().push_back(
moab, map_coords, one_prism));
CHKERR moab.write_file(
"prism_mesh.vtk",
"VTK",
"", &meshset, 1);
CHKERR moab.write_file(
"one_prism_mesh.vtk",
"VTK",
"", &one_prism_meshset,
1);
}
return 0;
}
postProc(post_proc), mapCoords(map_coords) {}
constexpr double def_val[] = {0, 0, 0};
case MBVERTEX:
case MBEDGE:
case MBTRI:
case MBQUAD:
case MBPRISM:
break;
default:
}
if (
type == MBTRI && (side != 3 && side != 4))
if (
type == MBQUAD && (side == 3 || side == 4))
for (
int dd = 0;
dd != nb_dofs; ++
dd)
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
else
}
&node_coords(0, 0));
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
constexpr
double eps = 1e-12;
"Inconsistency between node coords and integration "
"points");
}
"Inconsistency between node coords and integration "
"points");
}
MB_TAG_CREAT | MB_TAG_DENSE, def_val);
}
auto to_str = [](
auto i) {
return boost::lexical_cast<std::string>(
i); };
std::string tag_name_base =
"PrismType" + to_str(
type) +
"Side" + to_str(side);
std::cout << tag_name_base << endl;
for (size_t rr = 0; rr != trans_base.size1(); ++rr) {
auto tag_name = tag_name_base + "Base" + to_str(rr);
MB_TAG_CREAT | MB_TAG_DENSE, def_val);
&trans_base(rr, 0));
}
}
for (
int gg = 0; gg !=
triCoords.size1(); ++gg)
for (
int dd = 0;
dd != 2; ++
dd)
}
}
prism(prism) {}
int TriFE::getRule(
int order_row,
int order_col,
int order_data) {
return -1; }
int order_data) {
int side, sense, offset;
std::array<int, 2> swap = {0, 1};
if (side == 3)
swap = std::array<int, 2>{1, 0};
for (
int gg = 0; gg !=
triCoords.size1(); ++gg)
for (
int dd = 0;
dd != 2; ++
dd)
}
prism(prism) {}
int QuadFE::getRule(
int order_row,
int order_col,
int order_data) {
return -1; }
int order_data) {
int side, sense, offset;
constexpr double normal[3][2] = {{-1, 0}, {-0.5, 0.5}, {0, 1}};
constexpr double origin[3][2] = {{1, 0}, {1, 0}, {0, 0}};
constexpr int swap[3][2] = {{0, 1}, {0, 1}, {1, 0}};
int gg = 0;
for (size_t rr = 0; rr != edge_verts.size(); ++rr) {
const double edge_dist =
x * normal[side][0] + y * normal[side][1];
gaussPts(swap[side][0], gg) = edge_dist;
}
}
}
:
EdgeEle(m_field), edgeBlocks(edge_blocks),
prism(prism) {}
int EdgeFE::getRule(
int order_row,
int order_col,
int order_data) {
return -1; }
int order_data) {
int side, sense, offset;
if (side >= 3 && side <= 5) {
} else {
constexpr double normal[3][2] = {{1, 0}, {-0.5, 0.5}, {0, 1}};
constexpr double origin[3][2] = {{0, 1}, {1, 0}, {0, 0}};
constexpr int side_map[9] = {0, 1, 2, -1, -1, -1, 0, 1, 2};
int num_nodes;
&*coords.data().begin());
side = side_map[side];
gaussPts.resize(2, edge_verts.size(),
false);
for (size_t gg = 0; gg != edge_verts.size(); ++gg) {
const double edge_dist =
x * normal[side][0] + y * normal[side][1];
}
}
}
template <typename OP>
postProc(post_proc), mapCoords(map_coords), prism(prism) {}
template <typename OP>
constexpr double def_val[] = {0, 0, 0};
case MBVERTEX:
case MBEDGE:
case MBTRI:
case MBQUAD:
break;
default:
}
for (
int dd = 0;
dd != nb_dofs; ++
dd)
std::cerr <<
"Indices: " << data.
getIndices() << std::endl;
"Indicices/data inconsistency %3.1f != %d",
}
int side_prism, sense, offset;
CHKERR postProc.side_number(prism, fe_ent, side_prism, sense, offset);
} else
CHKERR postProc.side_number(prism, ent, side_prism, sense, offset);
auto &coords_at_pts = OP::getCoordsAtGaussPts();
const size_t nb_gauss_pts = coords_at_pts.size1();
nodeHandles.reserve(nb_gauss_pts);
nodeHandles.clear();
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
auto it = mapCoords.find(
t);
if (it != mapCoords.end())
nodeHandles.emplace_back(it->node);
else
}
}
auto to_str = [](
auto i) {
return boost::lexical_cast<std::string>(
i); };
std::string tag_name_base =
"FEType" + to_str(OP::getNumeredEntFiniteElementPtr()->getEntType()) +
"Type" + to_str(
type) +
"Side" + to_str(side_prism);
std::string tag_prism_name_base =
"PrismType" + to_str(
type) +
"Side" + to_str(side_prism);
MatrixDouble prism_base(trans_base.size1(), trans_base.size2());
if (trans_base.size2() != nodeHandles.size())
trans_base.size2(), nodeHandles.size());
for (size_t rr = 0; rr != trans_base.size1(); ++rr) {
std::string tag_name = tag_name_base + "Base";
int prism_mode_side;
CHKERR postProc.side_number(prism, node, prism_mode_side, sense, offset);
tag_name += to_str(prism_mode_side);
} else {
tag_name += to_str(rr);
}
std::cout << tag_name << endl;
CHKERR postProc.tag_get_handle(tag_name.c_str(), 1, MB_TYPE_DOUBLE,
th,
MB_TAG_CREAT | MB_TAG_DENSE, def_val);
CHKERR postProc.tag_set_data(
th, &nodeHandles[0], nodeHandles.size(),
&trans_base(rr, 0));
auto tag_prism_name = tag_prism_name_base + "Base" + to_str(rr);
Tag th_prism;
CHKERR postProc.tag_get_handle(tag_prism_name.c_str(), th_prism);
CHKERR postProc.tag_get_data(th_prism, &nodeHandles[0],
nodeHandles.size(), &prism_base(rr, 0));
}
}
double s = 0;
for (
unsigned int ii = 0; ii <
m.size1(); ii++) {
for (
unsigned int jj = 0; jj <
m.size2(); jj++) {
s += std::abs(
m(ii, jj));
}
}
return s;
};
prism_base -= trans_base;
constexpr
double eps = 1e-6;
cout << "Inconsistent base " << tag_prism_name_base << " "
<< tag_name_base << " sum " << sum << endl;
"Inconsistent base %s sum %6.4e", tag_prism_name_base.c_str(),
sum);
}
}