-
Christophe Geuzaine authoredChristophe Geuzaine authored
physical_element.cpp 8.94 KiB
#include <iostream>
#include <cassert>
#include "physical_element.h"
#include "gmsh_io.h"
#include "sgr.hpp"
physical_element::physical_element()
{}
size_t
physical_element::original_position(void) const
{
return m_original_position;
}
int
physical_element::dimension(void) const
{
return m_dimension;
}
int
physical_element::orientation(void) const
{
return m_orientation;
}
double
physical_element::measure(void) const
{
return m_measure;
}
vec_quadpt_3d
physical_element::integration_points(void) const
{
return m_phys_quadpoints;
}
int
physical_element::element_tag() const
{
return m_element_tag;
}
int
physical_element::element_type() const
{
return m_gmsh_elemtype;
}
vec_size_t
physical_element::node_tags(void) const
{
return m_node_tags;
}
/* Return only one determinant: this is used for geometric order 1
* where the determinants and the jacobians are constant */
double
physical_element::determinant(void) const
{
assert(m_geometric_order == 1);
return m_determinants[0];
}
double
physical_element::determinant(size_t iQp) const
{ assert(iQp < m_determinants.size());
return m_determinants[iQp];
}
/* Return all the determinants in all the integration points: this
* is for elements with geometric order > 1 */
vec_double
physical_element::determinants(void) const
{
assert(m_geometric_order >= 1);
return m_determinants;
}
/* Return only one jacobian: this is used for geometric order 1
* where the determinants and the jacobians are constant */
mat3d
physical_element::jacobian(void) const
{
assert(m_geometric_order == 1);
return m_jacobians[0];
}
/* Return the iQp-th jacobian: this is used for all
* geometric orders 1 */
mat3d
physical_element::jacobian(size_t iQp) const
{
assert(iQp < m_jacobians.size());
return m_jacobians[iQp];
}
/* Return all the jacobians in all the integration points: this
* is for elements with geometric order > 1 */
vec_mat3d
physical_element::jacobians(void) const
{
assert(m_geometric_order >= 1);
return m_jacobians;
}
point_3d
physical_element::barycenter(void) const
{
return m_barycenter;
}
physical_elements_factory::physical_elements_factory(const entity_params& ep)
: dim(ep.dim), tag(ep.tag), elemType(ep.etype), geom_order(ep.gorder),
approx_order(ep.aorder)
{
gmm::getElementsByType(elemType, cellTags, cellNodesTags, tag);
gmm::getIntegrationPoints(elemType, quadrature_name(2*approx_order), cellIps, cellIws);
gmm::getJacobians(elemType, cellIps, cellJacs, cellDets, cellPpts, tag);
gmm::getBasisFunctionsOrientationForElements(elemType, basis_func_name(approx_order),
cellOrientations, tag);
assert(cellOrientations.size() == cellTags.size());
gmm::getBarycenters(elemType, tag, false, false, barycenters);
std::vector<double> coord;
#ifdef USE_INITIAL_4_8_4_API
gmm::getKeysForElements(elemType, basis_func_name(approx_order),
keypairs, coord, tag, false);
#else
gmm::getKeysForElements(elemType, basis_func_name(approx_order),
tagKeys, entityKeys, coord, tag, false);
assert(tagKeys.size() == entityKeys.size());
#endif
keys_per_elem = gmm::getNumberOfKeysForElements(elemType, basis_func_name(approx_order));#ifdef USE_INITIAL_4_8_4_API
assert(keys_per_elem*cellTags.size() == keypairs.size());
#else
assert(keys_per_elem*cellTags.size() == tagKeys.size());
#endif
}
std::vector<physical_element>
physical_elements_factory::get_elements()
{
auto num_elems = cellTags.size();
auto num_nodes = cellNodesTags.size()/cellTags.size();
assert( cellTags.size()*num_nodes == cellNodesTags.size() );
auto num_gf = cellIws.size();
assert( num_gf * num_elems == cellDets.size() );
std::vector<physical_element> ret;
ret.reserve(num_elems);
for (size_t elem = 0; elem < num_elems; elem++)
{
physical_element new_pe;
new_pe.m_geometric_order = geom_order;
new_pe.m_approximation_order = approx_order;
new_pe.m_original_position = elem;
new_pe.m_dimension = dim;
new_pe.m_parent_entity_tag = tag;
new_pe.m_orientation = cellOrientations[elem];
new_pe.m_element_tag = cellTags[elem];
new_pe.m_measure = 0.0;
new_pe.m_gmsh_elemtype = elemType;
new_pe.m_barycenter = point_3d(barycenters[3*elem+0],
barycenters[3*elem+1],
barycenters[3*elem+2]);
new_pe.m_bf_keys.resize(keys_per_elem);
for (size_t i = 0; i < keys_per_elem; i++)
{
#ifdef USE_INITIAL_4_8_4_API
auto [vi, vu] = keypairs[keys_per_elem*elem + i];
#else
auto vi = tagKeys[keys_per_elem*elem + i];
auto vu = entityKeys[keys_per_elem*elem + i];
#endif
new_pe.m_bf_keys[i] = bf_key(vi,vu);
}
new_pe.m_node_tags.resize(num_nodes);
for (size_t i = 0; i < num_nodes; i++)
new_pe.m_node_tags[i] = cellNodesTags[num_nodes*elem+i];
auto elem_base = elem*num_gf;
for (size_t gf = 0; gf < num_gf; gf++)
{
auto gf_offset = elem_base + gf;
const auto JSIZE = 3*3; /* Jacobian matrix size */
auto jacs_base = gf_offset*JSIZE;
assert(jacs_base+8 < cellJacs.size());
/* GMSH returns jacobians by column */
Eigen::Matrix3d jac;
jac(0,0) = cellJacs[jacs_base+0]; /* dx/du */
jac(1,0) = cellJacs[jacs_base+1]; /* dy/du */
jac(2,0) = cellJacs[jacs_base+2]; /* dz/du */
jac(0,1) = cellJacs[jacs_base+3]; /* dx/dv */
jac(1,1) = cellJacs[jacs_base+4]; /* dy/dv */
jac(2,1) = cellJacs[jacs_base+5]; /* dz/dv */ jac(0,2) = cellJacs[jacs_base+6]; /* dx/dw */
jac(1,2) = cellJacs[jacs_base+7]; /* dy/dw */
jac(2,2) = cellJacs[jacs_base+8]; /* dz/dw */
new_pe.m_jacobians.push_back(jac);
auto dets_ofs = gf_offset;
assert( dets_ofs < cellDets.size() );
new_pe.m_determinants.push_back(cellDets[dets_ofs]);
new_pe.m_measure += cellDets[dets_ofs]*cellIws[gf];
const auto PSIZE = 3;
auto pts_base = gf_offset*PSIZE;
assert(pts_base + 2 < cellPpts.size());
point_3d pt(cellPpts[pts_base+0], cellPpts[pts_base+1], cellPpts[pts_base+2]);
quadrature_point_3d qpt(pt, cellIws[gf]*cellDets[dets_ofs]);
new_pe.m_phys_quadpoints.push_back(qpt);
}
ret.push_back( std::move(new_pe) );
}
return ret;
}
element_key::element_key()
{}
element_key::element_key(const physical_element& pe)
{
m_dim = pe.dimension();
m_elemType = pe.element_type();
auto nodes = pe.node_tags();
if (m_dim == 2)
{
m_key_data[0] = 3;
m_key_data[1] = nodes[0];
m_key_data[2] = nodes[1];
m_key_data[3] = nodes[2];
std::sort(std::next(m_key_data.begin()), m_key_data.end());
return;
}
std::stringstream ss;
ss << "Key not implemented for element type " << m_elemType;
throw std::invalid_argument(ss.str());
}
bool
element_key::operator<(const element_key& other) const
{
if (m_dim != other.m_dim or m_elemType != other.m_elemType)
return false;
if (m_key_data[0] != other.m_key_data[0])
return false;
auto mb = std::next(m_key_data.begin());
auto me = m_key_data.end();
auto ob = std::next(other.m_key_data.begin());
auto oe = other.m_key_data.end();
return std::lexicographical_compare(mb, me, ob, oe);
}
bool
element_key::operator==(const element_key& other) const
{
if (m_dim != other.m_dim or m_elemType != other.m_elemType) return false;
if (m_key_data[0] != other.m_key_data[0])
return false;
auto mb = std::next(m_key_data.begin());
auto me = m_key_data.end();
auto ob = std::next(other.m_key_data.begin());
return std::equal(mb, me, ob);
}
std::ostream& operator<<(std::ostream& os, const element_key& ek)
{
os << "ek: " << ek.m_dim << " " << ek.m_elemType << " (";
for (size_t i = 0; i < 8; i++)
os << ek.m_key_data[i] << " ";
os << ")";
return os;
}
element_key_factory::element_key_factory(int dim, int tag, int etype)
{
assert(dim == 2);
std::vector<size_t> nTags;
gmm::getElementFaceNodes(etype, 3, nTags, tag, true);
for (size_t i = 0; i < nTags.size(); i+= 3)
{
element_key ek;
ek.m_dim = dim;
ek.m_elemType = etype;
ek.m_key_data[0] = 3;
ek.m_key_data[1] = nTags[i+0];
ek.m_key_data[2] = nTags[i+1];
ek.m_key_data[3] = nTags[i+2];
std::sort(std::next(ek.m_key_data.begin()), ek.m_key_data.end());
ekeys.push_back( ek );
}
std::sort(ekeys.begin(), ekeys.end());
/* Not sure why there are duplicates */
ekeys.erase( std::unique(ekeys.begin(), ekeys.end()), ekeys.end() );
}
std::vector<element_key>
element_key_factory::get_keys(void) const
{
return ekeys;
}