Select Git revision
entity.cpp 23.79 KiB
#include <iostream>
#include <cassert>
#include "sgr.hpp"
#include "gmsh_io.h"
#include "entity.h"
#define NUM_TET_FACES 4
entity::entity(const entity_params& ep)
: dim(ep.dim), tag(ep.tag), elemType(ep.etype), g_order(ep.gorder),
a_order(ep.aorder), cur_elem_ordering(entity_ordering::GMSH),
m_dof_base(0), m_flux_base(0), m_index_base(0)
{
/* Prepare reference elements */
reference_elements_factory ref(ep);
reference_cells = ref.get_elements();
physical_elements_factory pef(ep);
physical_cells = pef.get_elements();
assert(dim == 3); /* otherwise you need getElementEdgeNodes() */
gmm::getElementFaceNodes(elemType, 3, faceNodesTags, tag); // 3 is for triangles
int tag2D = gm::addDiscreteEntity(dim-1);
int elemType2D = gmm::getElementType("triangle", g_order);
gmm::addElementsByType(tag2D, elemType2D, {}, faceNodesTags);
faceNodesTags.clear();
gmm::getElementsByType(elemType2D, faceTags, faceNodesTags, tag2D);
entity_params ep_2D;
ep_2D.dim = 2;
ep_2D.tag = tag2D;
ep_2D.etype = elemType2D;
ep_2D.gorder = g_order;
ep_2D.aorder = a_order;
reference_elements_factory ref_2D(ep_2D);
reference_faces = ref_2D.get_elements();
physical_elements_factory pef_2D(ep_2D);
physical_faces = pef_2D.get_elements();
assert(physical_faces.size() == 4*physical_cells.size());
assert(faceTags.size() == 4*physical_cells.size());
}
std::vector<size_t>
entity::face_tags() const
{
return faceTags;
}
double
entity::measure(void) const
{
auto meas = 0.0;
for (auto& pe : physical_cells)
meas += pe.measure();
return meas;
}
/* Project a function on this entity and return the
* corresponding vector of dofs on this entity. */
vecxd
entity::project(const scalar_function& function) const
{
auto num_bf = reference_cells[0].num_basis_functions();
vecxd ret = vecxd::Zero( physical_cells.size() * num_bf ); for (size_t iT = 0; iT < physical_cells.size(); iT++)
{
const auto& pe = physical_cells[iT];
//assert( pe.orientation() == 0 );
assert( pe.orientation() < reference_cells.size() );
const auto& re = reference_cells[pe.orientation()];
const auto pqps = pe.integration_points();
const auto dets = pe.determinants();
assert(pqps.size() == dets.size());
matxd mass = matxd::Zero(num_bf, num_bf);
vecxd rhs = vecxd::Zero(num_bf);
for (size_t iQp = 0; iQp < pqps.size(); iQp++)
{
const auto& pqp = pqps[iQp];
const vecxd phi = re.basis_functions(iQp);
mass += pqp.weight() * phi * phi.transpose();
rhs += pqp.weight() * function(pqp.point()) * phi;
}
ret.segment(iT*num_bf, num_bf) = mass.ldlt().solve(rhs);
}
return ret;
}
/* Project a function on this entity and store the partial result in
* the global vector pf at the offset corresponding to this entity. */
void
entity::project(const scalar_function& function, vecxd& pf) const
{
auto num_bf = reference_cells[0].num_basis_functions();
for (size_t iT = 0; iT < physical_cells.size(); iT++)
{
const auto& pe = physical_cells[iT];
assert( pe.orientation() < reference_cells.size() );
const auto& re = reference_cells[pe.orientation()];
const auto pqps = pe.integration_points();
const auto dets = pe.determinants();
assert(pqps.size() == dets.size());
matxd mass = matxd::Zero(num_bf, num_bf);
vecxd rhs = vecxd::Zero(num_bf);
for (size_t iQp = 0; iQp < pqps.size(); iQp++)
{
const auto& pqp = pqps[iQp];
const vecxd phi = re.basis_functions(iQp);
mass += pqp.weight() * phi * phi.transpose();
rhs += pqp.weight() * function(pqp.point()) * phi;
}
pf.segment(m_dof_base + iT*num_bf, num_bf) = mass.ldlt().solve(rhs);
}
}
/* Project a function on this entity and store the partial result in
* the global vector pf at the offset corresponding to this entity. */
void
entity::project(const vector_function& function, vecxd& pfx, vecxd& pfy, vecxd& pfz) const
{
auto num_bf = reference_cells[0].num_basis_functions();
for (size_t iT = 0; iT < physical_cells.size(); iT++)
{
const auto& pe = physical_cells[iT];
assert( pe.orientation() < reference_cells.size() );
const auto& re = reference_cells[pe.orientation()]; const auto pqps = pe.integration_points();
const auto dets = pe.determinants();
assert(pqps.size() == dets.size());
matxd mass = matxd::Zero(num_bf, num_bf);
vecxd rhs_x = vecxd::Zero(num_bf);
vecxd rhs_y = vecxd::Zero(num_bf);
vecxd rhs_z = vecxd::Zero(num_bf);
for (size_t iQp = 0; iQp < pqps.size(); iQp++)
{
const auto& pqp = pqps[iQp];
const vecxd phi = re.basis_functions(iQp);
mass += pqp.weight() * phi * phi.transpose();
vec3d fval = function(pqp.point());
rhs_x += pqp.weight() * fval(0) * phi;
rhs_y += pqp.weight() * fval(1) * phi;
rhs_z += pqp.weight() * fval(2) * phi;
}
Eigen::LDLT<matxd> mass_ldlt(mass);
pfx.segment(m_dof_base + iT*num_bf, num_bf) = mass_ldlt.solve(rhs_x);
pfy.segment(m_dof_base + iT*num_bf, num_bf) = mass_ldlt.solve(rhs_y);
pfz.segment(m_dof_base + iT*num_bf, num_bf) = mass_ldlt.solve(rhs_z);
}
}
const reference_element&
entity::cell_refelem(size_t iO) const
{
assert(iO < reference_cells.size());
return reference_cells[iO];
}
const reference_element&
entity::cell_refelem(const physical_element& pe) const
{
assert(pe.orientation() < reference_cells.size());
return reference_cells[pe.orientation()];
}
size_t
entity::cell_global_index(size_t iT) const
{
assert(iT < physical_cells.size());
return m_index_base + iT;
}
size_t
entity::cell_global_index_by_gmsh(size_t iT) const
{
assert(iT < physical_cells.size());
return m_index_base + physical_cells[iT].original_position();
}
size_t
entity::cell_local_dof_offset(size_t iT) const
{
assert(reference_cells.size() > 0);
return iT * reference_cells[0].num_basis_functions();
}
size_t
entity::cell_global_dof_offset(size_t iT) const
{
assert(reference_cells.size() > 0);
return m_dof_base + iT * reference_cells[0].num_basis_functions();
}
const reference_element&
entity::face_refelem(const physical_element& pe) const
{
assert(pe.orientation() < reference_faces.size());
return reference_faces[pe.orientation()];
}
size_t
entity::num_cells(void) const
{
return physical_cells.size();
}
size_t
entity::num_faces(void) const
{
return physical_faces.size();
}
size_t
entity::num_cell_orientations(void) const
{
return reference_cells.size();
}
std::vector<size_t>
entity::num_cells_per_orientation(void) const
{
std::vector<size_t> ret( num_cell_orientations(), 0 );
for (auto& pe : physical_cells)
{
assert( pe.orientation() < ret.size() );
ret[ pe.orientation() ]++;
}
return ret;
}
const physical_element&
entity::cell(size_t iT) const
{
assert(iT < physical_cells.size());
return physical_cells[iT];
}
const physical_element&
entity::face(size_t iF) const
{
assert(iF < physical_faces.size());
return physical_faces[iF];
}
size_t
entity::num_dofs(void) const
{
assert(reference_cells.size() > 0);
return physical_cells.size() * reference_cells[0].num_basis_functions();
}
size_t
entity::num_fluxes(void) const
{
assert(reference_faces.size() > 0);
return physical_faces.size() * reference_faces[0].num_basis_functions();
}
void
entity::base(size_t d_base, size_t f_base, size_t i_base)
{
m_dof_base = d_base; m_flux_base = f_base;
m_index_base = i_base;
}
size_t
entity::dof_base(void) const
{
return m_dof_base;
}
size_t
entity::flux_base(void) const
{
return m_flux_base;
}
size_t
entity::index_base(void) const
{
return m_index_base;
}
matxd
entity::mass_matrix(size_t iT) const
{
const auto& pe = physical_cells[iT];
assert( pe.orientation() < reference_cells.size() );
const auto& re = reference_cells[pe.orientation()];
const auto rqps = re.integration_points();
const auto dets = pe.determinants();
auto num_bf = re.num_basis_functions();
matxd mass = matxd::Zero(num_bf, num_bf);
for (size_t iQp = 0; iQp < rqps.size(); iQp++)
{
const auto& qp = rqps[iQp];
const vecxd phi = re.basis_functions(iQp);
mass += dets[iQp] * qp.weight() * phi * phi.transpose();
}
return mass;
}
void
entity::sort_by_orientation(void)
{
if (cur_elem_ordering == entity_ordering::BY_ORIENTATION)
return;
/* Sort by orientation and keep GMSH ordering within the same orientation */
auto comp = [](const physical_element& e1, const physical_element& e2) -> bool {
auto oe1 = e1.orientation();
auto oe2 = e2.orientation();
auto op1 = e1.original_position();
auto op2 = e2.original_position();
if ( oe1 < oe2 or (oe1 == oe2 and op1 < op2) )
return true;
return false;
};
/* Sort cells */
std::sort(physical_cells.begin(), physical_cells.end(), comp);
/* Sort all the stuff related to faces */
std::vector<physical_element> new_pf;
new_pf.resize( physical_faces.size() );
std::vector<size_t> new_ft; new_ft.resize( faceTags.size() );
std::vector<size_t> new_fnt;
new_fnt.resize( faceNodesTags.size() );
size_t nodes_per_face = faceNodesTags.size()/faceTags.size();
assert(nodes_per_face * faceTags.size() == faceNodesTags.size());
for (size_t new_iT = 0; new_iT < physical_cells.size(); new_iT++)
{
auto old_iT = physical_cells[new_iT].original_position();
for (size_t iF = 0; iF < NUM_TET_FACES; iF++)
{
auto old_ofs = NUM_TET_FACES*old_iT + iF;
auto new_ofs = NUM_TET_FACES*new_iT + iF;
new_pf[new_ofs] = physical_faces[old_ofs];
new_ft[new_ofs] = faceTags[old_ofs];
for (size_t iN = 0; iN < nodes_per_face; iN++)
{
auto old_ofs_nt = nodes_per_face*old_ofs + iN;
auto new_ofs_nt = nodes_per_face*new_ofs + iN;
new_fnt[new_ofs_nt] = faceNodesTags[old_ofs_nt];
}
}
}
std::swap(new_pf, physical_faces);
std::swap(new_ft, faceTags);
std::swap(new_fnt, faceNodesTags);
cur_elem_ordering = entity_ordering::BY_ORIENTATION;
}
void
entity::sort_by_gmsh(void)
{
/* Sort by GMSH original ordering */
auto comp = [](const physical_element& e1, const physical_element& e2) -> bool {
auto op1 = e1.original_position();
auto op2 = e2.original_position();
if ( op1 < op2 )
return true;
return false;
};
std::sort(physical_cells.begin(), physical_cells.end(), comp);
cur_elem_ordering = entity_ordering::GMSH;
throw std::logic_error("NOT SORTED");
}
int
entity::gmsh_tag(void) const
{
return tag;
}
int
entity::gmsh_elem_type(void) const
{
return elemType;
}
entity_ordering
entity::current_elem_ordering(void) const
{
return cur_elem_ordering;
}
void
entity::populate_differentiation_matrices_planar(entity_data_cpu& ed) const
{
/* In the case of geometric order = 1, the inverse of the
* mass matrix is embedded in the differentiation matrix */
assert(g_order == 1 and g_order == ed.g_order);
auto dm_rows = ed.num_orientations * ed.num_bf;
auto dm_cols = 3 * ed.num_bf;
ed.differentiation_matrices = matxd::Zero(dm_rows, dm_cols);
for (size_t iO = 0; iO < ed.num_orientations; iO++)
{
matxd dm = matxd::Zero(ed.num_bf, 3*ed.num_bf);
auto& re = reference_cells[iO];
auto qps = re.integration_points();
auto num_qp = qps.size();
for (size_t iQp = 0; iQp < num_qp; iQp++)
{
auto w = qps[iQp].weight();
vecxd phi = re.basis_functions(iQp);
matxd dphi = re.basis_gradients(iQp);
vecxd du_phi = dphi.col(0);
vecxd dv_phi = dphi.col(1);
vecxd dw_phi = dphi.col(2);
dm.block(0, 0*ed.num_bf, ed.num_bf, ed.num_bf) += w * phi * du_phi.transpose();
dm.block(0, 1*ed.num_bf, ed.num_bf, ed.num_bf) += w * phi * dv_phi.transpose();
dm.block(0, 2*ed.num_bf, ed.num_bf, ed.num_bf) += w * phi * dw_phi.transpose();
}
matxd mass = re.mass_matrix();
Eigen::LDLT<matxd> mass_ldlt(mass);
if (mass_ldlt.info() != Eigen::Success )
throw std::invalid_argument("LDLT failed");
size_t base_row = iO * ed.num_bf;
ed.differentiation_matrices.block(base_row, 0, ed.num_bf, 3*ed.num_bf) =
mass_ldlt.solve(dm);
}
}
void
entity::populate_differentiation_matrices_curved(entity_data_cpu& ed) const
{
assert(g_order > 1 and g_order == ed.g_order);
auto dm_rows = ed.num_orientations * ed.num_bf * ed.num_qp;
auto dm_cols = 3 * ed.num_bf;
ed.differentiation_matrices = matxd::Zero(dm_rows, dm_cols);
for (size_t iO = 0; iO < ed.num_orientations; iO++)
{
auto& re = reference_cells[iO];
auto qps = re.integration_points();
auto num_qp = qps.size();
for (size_t iQp = 0; iQp < num_qp; iQp++)
{
auto w = qps[iQp].weight();
vecxd phi = re.basis_functions(iQp);
matxd dphi = re.basis_gradients(iQp);
vecxd du_phi = dphi.col(0);
vecxd dv_phi = dphi.col(1);
vecxd dw_phi = dphi.col(2);
auto row_ofs = ed.num_bf*ed.num_qp*iO + ed.num_bf*iQp;
ed.differentiation_matrices.block(row_ofs, 0*ed.num_bf, ed.num_bf, ed.num_bf) =
w * phi * du_phi.transpose();
ed.differentiation_matrices.block(row_ofs, 1*ed.num_bf, ed.num_bf, ed.num_bf) =
w * phi * dv_phi.transpose();
ed.differentiation_matrices.block(row_ofs, 2*ed.num_bf, ed.num_bf, ed.num_bf) =
w * phi * dw_phi.transpose();
}
}
}
void
entity::populate_differentiation_matrices(entity_data_cpu& ed) const
{
assert(g_order == ed.g_order);
if (ed.g_order == 1)
populate_differentiation_matrices_planar(ed);
else
populate_differentiation_matrices_curved(ed);
}
static size_t face_orientations[] = {
1, 0, 1, 4, 1, 0, 0, 2, 0, 0, 1, 5,
1, 1, 0, 0, 0, 1, 1, 3, 0, 1, 0, 1,
3, 2, 1, 4, 3, 2, 0, 2, 2, 0, 3, 5,
1, 3, 2, 0, 2, 1, 3, 3, 0, 3, 2, 1,
5, 2, 3, 4, 3, 4, 2, 2, 4, 2, 3, 5,
3, 5, 2, 0, 2, 3, 5, 3, 2, 3, 4, 1,
5, 4, 5, 4, 5, 4, 4, 2, 4, 4, 5, 5,
5, 5, 4, 0, 4, 5, 5, 3, 4, 5, 4, 1
};
#define FACE_ORIENTATION(co, fn) ( face_orientations[4*co + fn] )
void
entity::populate_lifting_matrices_planar(entity_data_cpu& ed,
const dofs_f2c_t& dofs_f2c) const
{
assert(g_order == 1 and g_order == ed.g_order);
auto lm_rows = ed.num_orientations * ed.num_bf;
auto lm_cols = NUM_TET_FACES * ed.num_fluxes;
ed.lifting_matrices = matxd::Zero(lm_rows, lm_cols);
for (size_t iO = 0; iO < ed.num_orientations; iO++)
{
if (dofs_f2c.find(iO) == dofs_f2c.end())
continue;
auto& f2c = dofs_f2c.at(iO);
assert(f2c.size() == NUM_TET_FACES*ed.num_fluxes);
matxd lm = matxd::Zero(ed.num_bf, lm_cols);
for (size_t iF = 0; iF < NUM_TET_FACES; iF++)
{
auto fO = FACE_ORIENTATION(iO, iF);
assert(fO < reference_faces.size());
auto &rf = reference_faces[fO];
matxd mm = rf.mass_matrix();
assert(mm.rows() == ed.num_fluxes);
for (size_t i = 0; i < mm.rows(); i++)
{
assert(iF*ed.num_fluxes + i < f2c.size());
auto lm_row = f2c[iF*ed.num_fluxes + i];
assert(lm_row < ed.num_bf);
for (size_t j = 0; j < mm.cols(); j++)
{
auto lm_col = iF*ed.num_fluxes + j;
assert(lm_col < lm_cols);
lm(lm_row, lm_col) = mm(i, j); }
}
}
auto& re = reference_cells[iO];
matxd mass = re.mass_matrix();
Eigen::LDLT<matxd> mass_ldlt(mass);
if (mass_ldlt.info() != Eigen::Success )
throw std::invalid_argument("LDLT failed");
size_t base_row = iO * ed.num_bf;
ed.lifting_matrices.block(base_row, 0, ed.num_bf, lm_cols) =
mass_ldlt.solve(lm);
//std::cout << "Orientation : " << iO << std::endl;
//std::cout << mass_ldlt.solve(lm) << std::endl;
}
}
void
entity::populate_lifting_matrices_curved(entity_data_cpu&, const dofs_f2c_t&) const
{
}
void
entity::populate_lifting_matrices(entity_data_cpu& ed) const
{
assert(reference_faces.size() > 0);
auto num_face_bf = reference_faces[0].num_basis_functions();
dofs_f2c_t dofs_f2c;
for (size_t iT = 0; iT < physical_cells.size(); iT++)
{
auto& pe = physical_cells[iT];
auto orient = pe.orientation();
auto cnt = pe.bf_keys();
#ifndef EXPENSIVE_ASSERTS
if ( dofs_f2c.find(orient) != dofs_f2c.end() )
continue;
#endif
std::vector<size_t> f2c;
f2c.reserve( NUM_TET_FACES*num_face_bf );
for (size_t iF = 0; iF < NUM_TET_FACES; iF++)
{
auto face_ofs = NUM_TET_FACES*iT+iF;
auto& pf = physical_faces[face_ofs];
auto fnt = pf.bf_keys();
for (size_t i = 0; i < fnt.size(); i++)
for (size_t j = 0; j < cnt.size(); j++)
if (fnt[i] == cnt[j])
{
assert(f2c.size() == iF*num_face_bf + i);
f2c.push_back(j);
}
}
assert(f2c.size() == NUM_TET_FACES*num_face_bf);
#ifdef EXPENSIVE_ASSERTS
if ( dofs_f2c.find(orient) != dofs_f2c.end() )
assert(dofs_f2c[orient] == f2c);
#endif
dofs_f2c[orient] = f2c;
}
if (g_order == 1)
populate_lifting_matrices_planar(ed, dofs_f2c);
else
populate_lifting_matrices_curved(ed, dofs_f2c);
}
void
entity::populate_jacobians(entity_data_cpu& ed) const
{
auto num_elems = num_cells();
if (g_order == 1)
{
ed.jacobians = matxd::Zero(3*num_elems, 3);
ed.cell_determinants = vecxd::Zero(num_elems);
for (size_t iT = 0; iT < num_elems; iT++)
{
auto& pe = physical_cells[iT];
ed.jacobians.block(3*iT, 0, 3, 3) = pe.jacobian().inverse().transpose();
ed.cell_determinants(iT) = pe.determinant();
}
assert(num_faces() == 4*num_elems);
ed.face_determinants = vecxd(4*num_elems);
for (size_t iF = 0; iF < num_faces(); iF++)
{
auto& pf = physical_faces[iF];
ed.face_determinants(iF) = pf.determinant();
}
}
else
{
ed.jacobians = matxd::Zero(3*num_elems*ed.num_qp, 3);
ed.cell_determinants = vecxd::Zero(num_elems*ed.num_qp);
for (size_t iT = 0; iT < num_elems; iT++)
{
for (size_t iQp = 0; iQp < ed.num_qp; iQp++)
{
auto& pe = physical_cells[iT];
ed.jacobians.block(3*iT*ed.num_qp + 3*iQp, 0, 3, 3) =
pe.jacobian(iQp).inverse().transpose();
ed.cell_determinants(iT*ed.num_qp + iQp) = pe.determinant(iQp);
}
}
}
}
void
entity::populate_normals(entity_data_cpu& ed) const
{
size_t num_face_qps = 1;
size_t num_elems = num_cells();
if (g_order == 1)
ed.normals = matxd::Zero(NUM_TET_FACES*num_elems, 3);
else
{
assert(reference_faces.size() > 0);
num_face_qps = reference_faces[0].num_integration_points();
ed.normals = matxd::Zero(NUM_TET_FACES*num_elems*ed.num_qp, 3);
}
Eigen::Vector3d ref_normal;
ref_normal(0) = 0; ref_normal(1) = 0; ref_normal(2) = 1;
for (size_t iT = 0; iT < num_cells(); iT++)
{
for (size_t iF = 0; iF < NUM_TET_FACES; iF++)
{
for (size_t iQp = 0; iQp < num_face_qps; iQp++)
{ auto fnum = NUM_TET_FACES*iT+iF;
assert(fnum < physical_faces.size());
Eigen::Vector3d n = physical_faces[fnum].jacobian(iQp)*ref_normal;
auto norm_offset = (NUM_TET_FACES*iT+iF)*num_face_qps + iQp;
assert(norm_offset < ed.normals.rows());
ed.normals.block(norm_offset, 0, 1, 3) = n.transpose()/n.norm();
}
}
}
}
void
entity::populate_entity_data(entity_data_cpu& ed) const
{
assert(cur_elem_ordering == entity_ordering::BY_ORIENTATION);
assert(reference_faces.size() > 0);
assert(reference_cells.size() > 0);
ed.ordering = cur_elem_ordering;
ed.g_order = g_order;
ed.a_order = a_order;
ed.num_elems = num_cells_per_orientation();
ed.num_orientations = num_cell_orientations();
ed.num_bf = reference_cells[0].num_basis_functions();
ed.num_fluxes = reference_faces[0].num_basis_functions();
ed.num_qp = (g_order == 1) ? 1 : reference_cells[0].num_integration_points();
ed.num_face_qp = (g_order == 1) ? 1 : reference_faces[0].num_integration_points();
ed.dof_base = m_dof_base;
ed.flux_base = m_flux_base;
populate_differentiation_matrices(ed);
populate_jacobians(ed);
populate_normals(ed);
populate_lifting_matrices(ed);
ed.invmass_matrices = matxd::Zero(ed.num_bf*num_cells(), ed.num_bf);
for (size_t iT = 0; iT < num_cells(); iT++)
{
matxd mass = mass_matrix(iT);
ed.invmass_matrices.block(iT*ed.num_bf, 0, ed.num_bf, ed.num_bf) =
mass.inverse();
}
ed.fluxdofs_mine.resize(4*ed.num_fluxes*num_cells());
for (size_t iT = 0; iT < num_cells(); iT++)
{
auto& pe = physical_cells[iT];
auto cell_keys = pe.bf_keys();
auto cell_base = iT * ed.num_bf;
for (size_t iF = 0; iF < 4; iF++)
{
auto face_num = 4*iT+iF;
auto& pf = physical_faces[face_num];
auto face_keys = pf.bf_keys();
auto face_base = face_num * ed.num_fluxes;
for (size_t i = 0; i < face_keys.size(); i++)
for (size_t j = 0; j < cell_keys.size(); j++)
if ( face_keys.at(i) == cell_keys.at(j) )
ed.fluxdofs_mine.at(face_base + i) = cell_base + j;
}
}
}
std::pair<std::vector<point_3d>, std::vector<quadrature_point_3d>>
integrate(const entity& e, int order)
{
auto tag = e.gmsh_tag();
auto elemType = e.gmsh_elem_type();
std::vector<double> ips;
std::vector<double> iws;
gmm::getIntegrationPoints(elemType, quadrature_name(order), ips, iws);
auto num_qp = iws.size();
std::vector<double> ppts;
std::vector<double> jacs;
std::vector<double> dets;
gmm::getJacobians(elemType, ips, jacs, dets, ppts, tag);
std::vector<point_3d> retr;
for (size_t iQp = 0; iQp < num_qp; iQp++)
{
double u = ips[3*iQp + 0];
double v = ips[3*iQp + 1];
double w = ips[3*iQp + 2];
point_3d p(u,v,w);
retr.push_back(p);
}
std::vector<quadrature_point_3d> retp;
retp.reserve(dets.size());
for (auto& pe : e)
{
auto g_pos = pe.original_position();
auto dets_base = g_pos*num_qp;
auto ppts_base = 3*g_pos*num_qp;
for (size_t iQp = 0; iQp < num_qp; iQp++)
{
double w = dets[dets_base+iQp]*iws[iQp];
double x = ppts[ppts_base + 3*iQp + 0];
double y = ppts[ppts_base + 3*iQp + 1];
double z = ppts[ppts_base + 3*iQp + 2];
point_3d p(x,y,z);
quadrature_point_3d q(p,w);
retp.push_back(q);
}
}
return std::make_pair(retr, retp);
}