Homec4science
Diffusion CADD_mesher (dev/incorporate-qhull) Browse

rCADDMESH/src/domain/cadd_mesh.cc
58a104d7bd95dev/incorporate-qhull

Last Change
Raw File

cadd_mesh.cc
View Options

extern "C" {
#include <stdio.h>
#include <libqhull/libqhull.h>
#include <libqhull/mem.h>
#include <libqhull/qset.h>
#include <libqhull/poly.h> /* for qh_vertexneighbors() */
#include <libqhull/geom.h> /* for qh_facetarea(facet)*/
}
#include "cadd_mesh.hh"
#include <cmath>
#include <sstream>
#include <limits>
#include <set>
#include <algorithm>
#include <stdexcept>
#include <gsl/gsl_cblas.h>
#include "dumper_paraview.hh"
// akantu includes for mesh io
#include <aka_common.hh>
#include <aka_array.hh>
#include <mesh_io_msh.hh>
#include <mesh_utils.hh>
#include "mesh_helper.hh"
#include "mini_mesh.hh"
/* -------------------------------------------------------------------------- */
extern "C" {
// LU decomoposition of a general matrix
void dgetrf_(int* M, int *N, double* A, int* lda, int* IPIV, int* INFO);
// generate inverse of a matrix given its LU decomposition
void dgetri_(int* N, double* A, int* lda, int* IPIV, double* WORK, int* lwork, int* INFO);
}
/* -------------------------------------------------------------------------- */
void inverse(double * A, int N)
{
int IPIV[N+1];
int LWORK = N*N;
double WORK[LWORK];
int INFO;
dgetrf_(&N,&N,A,&N,IPIV,&INFO);
dgetri_(&N,A,&N,IPIV,WORK,&LWORK,&INFO);
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
Mesh<DIM>::Mesh(const std::vector<Real> & points, Real tol,
Geometry<DIM> * excluded_points)
:permutation(NULL), remutation(NULL), radius_ratio_computed(false),
meshed(false), cleaned(false), interpolation_initialised(false),
aka_mesh(NULL), aka_mesh_exists(false), tol(tol)
{
for (Uint i = 0; i < points.size() ; i += DIM) {
if ((excluded_points == NULL) ||
(!excluded_points->is_inside(&points.at(i),
this->tol*0)) ) {
this->renumerotation.push_back(i/DIM);
PointRef<DIM> tmp;
for (Uint j = 0; j < DIM ; ++j) {
tmp[j] = points.at(i+j);
}
this->nodes.addPoint(tmp);
}
}
std::stringstream number;
number << this->aka_mesh_counter++;
this->my_aka_mesh_name = this->aka_mesh_name.substr(0, 5) + number.str();
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
Mesh<DIM>::~Mesh() {
if (this->permutation != NULL) {
delete this->permutation;
delete this->remutation;
}
if (this->aka_mesh != NULL) {
delete this->aka_mesh;
}
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Mesh<DIM>::printself(std::ostream & stream, int indent) const
{
indent_space(stream, indent);
stream << "Mesh:" << std::endl;
if (this->meshed) {
indent_space(stream, indent+indent_incr);
stream << "has been meshed" << std::endl;
print_vector(this->nodes.getVector(), "Nodes", DIM, indent+indent_incr, stream);
print_vector(this->connectivity, "Connectivity", DIM+1, indent + indent_incr, stream);
} else {
indent_space(stream, indent+indent_incr);
stream << "has NOT been meshed" << std::endl;
}
if (this->cleaned) {
indent_space(stream, indent+indent_incr);
stream << "has been cleaned" << std::endl;
} else {
indent_space(stream, indent+indent_incr);
stream << "has NOT been cleaned" << std::endl;
}
if (this->interpolation_initialised) {
indent_space(stream, indent+indent_incr);
stream << "Is ready to be used for interpolations" << std::endl;
} else {
indent_space(stream, indent+indent_incr);
stream << "Is NOT ready for interpolations" << std::endl;
}
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
inline bool degenerate_finder(const std::vector<double> & points,
const std::vector<int> & element,
double volume,
double criterion)
{
Uint nb_nodes = points.size()/DIM;
const Real * points_ptr = &points.front();
const int * element_ptr = &element.front();
Real ratio = computeRadiusRatio<DIM>(points_ptr, element_ptr, nb_nodes);
// smaller ratio than criterion means element is good
return (ratio < criterion) && (ratio > 0);
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
inline void center_of_mass(const std::vector<double> & points,
const std::vector<int> & element,
PointRef<DIM> & center_mass) {
for (Uint i = 0; i < DIM; ++i) {
center_mass[i] = 0.;
}
for (Uint nd = 0; nd < DIM + 1; ++nd) {
for (Uint dir = 0; dir < DIM; ++dir) {
center_mass[dir] += points[element[nd] * DIM + dir];
}
}
for (Uint i = 0; i < DIM; ++i) {
center_mass[i] /= DIM+1;
}
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
inline void shear_points(std::vector<Real> & new_coords,
const std::vector<Real> & ref_coords,
Uint modified_dir,
Uint dir2,
Real amount) {
Uint nb_pts = new_coords.size()/DIM;
for (Uint id = 0; id < nb_pts ; ++id ) {
new_coords[DIM*id+modified_dir] +=
amount * ref_coords[DIM*id + dir2];
}
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
inline void smart_joggle(std::vector<Real> & new_coords,
const std::vector<Real> & ref_coords) {
for (Uint shear_dim = 0; shear_dim < DIM - 1; ++shear_dim) {
for (Uint ref_dim = shear_dim + 1; ref_dim < DIM; ++ref_dim) {
shear_points<DIM>(new_coords, ref_coords, shear_dim, ref_dim,
(ref_dim + shear_dim) * 1e-1);
}
}
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Mesh<DIM>::delaunay(const Real & degeneration_criterion,
const std::string & flags){
// make a copy of the nodes before they get jiggled
std::vector<Real> copy_nodes = this->nodes.getVector();
smart_joggle<DIM>(copy_nodes, this->nodes.getVector());
boolT ismalloc = False; /* True if qhull should free points in
qh_freeqhull() or reallocation */
int exitcode; /* 0 if no error from qhull */
//char flags[] = "qhull d QJ Pg Ft"; /* flags for qhull */
const Uint flagsize = 129;
char cflags[flagsize] = {"\0"};
if (flags.size() > flagsize - 1) {
std::stringstream error_stream ;
error_stream << "The flags you specified: '" << flags << "' are longer than the maximum size of " << flagsize - 1 << ". You probably fucked up and should reconsider your qhull options";
throw(std::out_of_range(error_stream.str()));
}
strcpy(cflags, flags.c_str()); /* flags for qhull */
/* Call qhull */
exitcode = qh_new_qhull(DIM, copy_nodes.size()/DIM, &copy_nodes[0],
ismalloc, cflags, NULL, stderr);
if (exitcode){
std::stringstream exit_stream;
exit_stream << "qhull failed with exitcode " << exitcode;
throw exit_stream.str();
} else {
PointContainer<DIM> dropped_points;
//I use an stl vector to build the connectivity list. making sure it's empty:
connectivity.clear();
//build neighbours:
qh_vertexneighbors();
//elements in 3d are 4d facets, hence
facetT *facet;
vertexT * vertex, ** vertexp;
Real element_volume;
int nb_elems = 0;
std::vector<int> current_element;
FORALLfacets {
if (facet->upperdelaunay)
continue;
element_volume = qh_facetarea(facet);
++nb_elems;
unsigned int counter = 0;
current_element.clear();
FOREACHvertex_(facet->vertices){
current_element.push_back(qh_pointid(vertex->point));
++ counter;
}
if (counter != DIM+1){
std::stringstream exit_stream;
exit_stream << "With flags '" << cflags << "': Element # " << nb_elems
<< " has " << counter << " vertices instead of "
<< DIM+1 << ". The vertices are ";
for (Uint vert = 0; vert < counter - 1; ++vert) {
exit_stream << current_element[vert] << ", ";
}
exit_stream << current_element.back() << "." << std::endl;
exit_stream << "Indices from:" << std::endl << this->nodes << std::endl;
throw exit_stream.str();
} else {
bool is_valid = degenerate_finder<DIM>(this->nodes.getVector(), current_element,
element_volume, degeneration_criterion);
if (is_valid) {
for (unsigned int i = 0; i < DIM+1; ++i) {
connectivity.push_back(current_element[i]);
}
} else {
PointRef<DIM> center;
center_of_mass(this->nodes.getVector(), current_element, center);
dropped_points.addPoint(center);
}
}
}
if (dropped_points.getSize() != 0) {
std::cerr << "WARNING: Elements have been dropped because they have been "
<< "evaluated and judged degenerated. Check the following list "
<< "of their centers of gravity carefully to make sure nothing "
<< "got fucked up!" << std::endl;
std::cerr << "Dropped elements with centers of gravity at " << dropped_points << std::endl;
}
}
/* Free the memory */
qh_freeqhull(qh_ALL);
this->meshed = true;
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Mesh<DIM>::cleanup(const Geometry<DIM> & atomistic_domain) {
Uint nb_points = this->nodes.getSize();
Uint nb_elems = this->connectivity.size()/(DIM+1);
std::vector<Real> new_points;
std::vector<int> new_connectivity;
this->permutation = new std::vector<int>(nb_points, -1);
this->remutation = new std::vector<int>(nb_points, -1);
Uint counter = 0;
// loop through elements
for (Uint i = 0; i < nb_elems ; ++i) {
Real center_gravity[DIM] = {0};
int * element_nodes;
try {
element_nodes = &connectivity.at(i*(DIM+1));
} catch (std::out_of_range & e) {
std::cout << "found the sucker. Tried to access element " << i*(DIM+1) << " but connectivity contains only " << connectivity.size() << "entries" << std::endl;
throw("scheisse");
}
Uint nb_pure_atoms = 0; // an element is invalid if all its nodes are free atoms
// loop through nodes of current element
for (Uint j = 0; j < DIM+1 ; ++j) {
// check for each node if it is valid
int point_id = element_nodes[j];
for (Uint k = 0; k < DIM ; ++k) {
center_gravity[k] += this->nodes[point_id][k];
}
if (atomistic_domain.is_inside(&this->nodes.getVector().at(DIM*point_id), 0.) ) {
nb_pure_atoms ++;
}
}
// if at least one node is not a pure atom, we're good to go
// also, if an element is valid, its nodes are (obviously) valid as well
if (nb_pure_atoms < DIM+1) {
// check the center of gravity
for (Uint k = 0; k < DIM ; ++k) {
center_gravity[k] /= (DIM+1);
}
if (!atomistic_domain.is_inside(center_gravity, 0.)) {
for (Uint elem_node = 0; elem_node < DIM+1 ; ++elem_node) {
int point_id = element_nodes[elem_node];
if (this->permutation->at(point_id) < 0) {
this->permutation->at(point_id) = counter;
this->remutation ->at(counter) = point_id;
++counter;
for (Uint k = 0; k < DIM ; ++k) {
new_points.push_back(this->nodes[point_id][k]);
}
}
new_connectivity.push_back(this->permutation->at(point_id));
}
}
}
}
this->nodes = PointContainer<DIM>(new_points);
this->connectivity = new_connectivity;
this->cleaned = true;
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Mesh<DIM>::initInterpolation() {
Uint nb_elems = this->connectivity.size()/(DIM+1);
for (Uint elem_id = 0; elem_id < nb_elems ; ++elem_id) {
int * element_nodes = &(this->connectivity.at(elem_id * (DIM + 1)));
// fill the first line of the transformation matrix is trivial, because only
// the first is 1
this->reverse_matrices.push_back(1.);
for (Uint i = 0; i < DIM ; ++i) {
this->reverse_matrices.push_back(0.);
}
// fill the matrix row by row
for (Uint dir = 0; dir < DIM ; ++dir) {
//Uint row = (dir + 1) * (DIM + 1);
Real last_node_coord = this->nodes[element_nodes[DIM]][dir];
this->reverse_matrices.push_back(last_node_coord);
for (Uint node_local = 0; node_local < DIM ; ++node_local) {
this->reverse_matrices.push_back(this->nodes[element_nodes[node_local]][dir]
- last_node_coord);
}
}
//invert the matrix to make it useful
inverse(&(this->reverse_matrices.at(elem_id*(DIM+1)*(DIM+1))), DIM + 1);
}
this->interpolation_initialised = true;
}
/* -------------------------------------------------------------------------- */
template<Uint DIM>
Real Mesh<DIM>::interpolate(const Real * coords,
const std::vector<Real> & nodal_values) const {
Real tol = 1e-6;
Uint nb_elems = this->connectivity.size()/(DIM+1);
Real elem_coords[DIM+2]={0}; //we'll use the last slot for the n-th shape function value
Real * shape_vals = elem_coords+1;
Real glob_coords[DIM+1];
glob_coords[0] = 1;
for (Uint dir = 0; dir < DIM ; ++dir) {
glob_coords[dir+1] = coords[dir];
}
Real max_viol = -std::numeric_limits<Real>::max();
Uint elem_id;
bool found = false;
int *element_nodes = NULL;
//loop through elements and figure out in which one the coords are
for (elem_id = 0; elem_id < nb_elems ; ++elem_id) {
const Real* matrix = &(this->reverse_matrices.at(elem_id * (DIM+1)*(DIM+1)));
cblas_dgemv(CblasRowMajor, CblasNoTrans, DIM+1, DIM+1, 1., matrix, DIM+1,
glob_coords, 1, 0., elem_coords, 1);
if (fabs(1-elem_coords[0]) > tol){
std::stringstream error_stream ;
error_stream << " Reversion matrix for element "<< elem_id
<< " is obviously wrong, because the sum of "
<< " Shape functions is not 1 " << std::endl;
throw(error_stream.str());
}
// check if any of the shape functions are < 0
Real minimum = std::numeric_limits<Real>::max();
Real sum = 0;
for (Uint xi_eta = 1; xi_eta < DIM+1 ; ++xi_eta) {
keep_min(minimum, elem_coords[xi_eta]);
sum += elem_coords[xi_eta];
}
shape_vals[DIM] = 1-sum;
keep_min(minimum, shape_vals[DIM]);
keep_max(max_viol, minimum);
if (minimum + tol >= 0) {
found = true;
element_nodes = const_cast<int *>(&(this->connectivity.at(elem_id*(DIM+1))));
break;
}
}
if (!found) {
std::stringstream error_stream;
error_stream << "It seems that the point (";
for (Uint i = 0; i < DIM-1 ; ++i) {
error_stream << coords[i] << ", ";
}
error_stream << coords[DIM-1] << ") is not within the domain. "
<< "the minimum constraint violation was " << max_viol
<< std::endl;
error_stream << "epsilon = " << tol
<< std::endl;
error_stream << "Is it possible that the auxiliary mesh is too small (e.g. "
<< "you specified its size in lattice constants instead of "
<< "distance units" << std::endl;
throw(error_stream.str());
}
Real interpolated_value = 0;
for (Uint i = 0; i < DIM+1 ; ++i) {
int index;
if (this->permutation != NULL) {
index = this->remutation->at(element_nodes[i]);
} else {
index = element_nodes[i];
}
try {
index = this->renumerotation.at(index);
} catch (std::out_of_range & e) {
throw (std::out_of_range("Trying to access a point index out of the range orf points. did you use the Qhull option QZ?"));
}
interpolated_value += nodal_values.at(index) * shape_vals[i];
}
return interpolated_value;
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Mesh<DIM>::dump_paraview(std::string filename) throw (std::string){
if (this->meshed) {
iohelper::DumperParaview dumper;
dumper.setMode(iohelper::TEXT);
// a copy on the altar of const correctness
std::vector<Real> tmp_nodes = this->nodes.getVector();
dumper.setPoints(&tmp_nodes[0], DIM, this->nodes.getSize(), filename);
iohelper::ElemType type;
if (DIM == 2) {
type = iohelper::TRIANGLE1;
} else if (DIM == 3) {
type = iohelper::TETRA1;
} else {
throw(std::string("Only 2 and 3D meshes can be dumped"));
}
dumper.setConnectivity(&this->connectivity[0], type, connectivity.size()/(DIM+1), iohelper::C_MODE);
dumper.setPrefix("./");
dumper.init();
dumper.dump();
} else {
throw(std::string("this mesh has not been delaunay'ed yet"));
}
}
/* -------------------------------------------------------------------------- */
inline void crossProduct (const Real * vec1, const Real * vec2,
Real * result){
for (Uint dir = 0 ; dir < threeD ; ++dir) {
Uint dir2 = (dir+1)%threeD;
Uint dir3 = (dir+2)%threeD;
result[dir] = vec1[dir2]*vec2[dir3] - vec1[dir3]*vec2[dir2];
}
}
/* -------------------------------------------------------------------------- */
template<Uint DIM>
Real signOfJacobian(const Uint * elem, const std::vector<Real> & nodes) {
if (DIM == 2) {
Real vec1[DIM], vec2[DIM];
for (Uint i = 0 ; i < DIM ; ++i) {
vec1[i] = nodes[DIM * elem[1] + i] - nodes [DIM * elem[0] +i];
vec2[i] = nodes[DIM * elem[2] + i] - nodes [DIM * elem[0] +i];
}
return vec1[0]*vec2[1] - vec1[1]*vec2[0];
} if (DIM == 3) {
Real vec1[DIM], vec2[DIM], vec3[DIM], tmp[DIM];
for (Uint i = 0 ; i < DIM ; ++i) {
vec1[i] = nodes[DIM * elem[1] + i] - nodes [DIM * elem[0] +i];
vec2[i] = nodes[DIM * elem[2] + i] - nodes [DIM * elem[0] +i];
vec3[i] = nodes[DIM * elem[3] + i] - nodes [DIM * elem[0] +i];
}
crossProduct(vec1, vec2, tmp);
// compute scalar product
Real val = 0;
for (Uint i = 0 ; i < DIM ; ++i) {
val += vec3[i] * tmp[i];
}
return val;
}
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Mesh<DIM>::create_aka_mesh() {
this->aka_mesh = new akantu::Mesh(DIM, this->my_aka_mesh_name, 0);
// fill the nodes of the aka_mesh
akantu::Array<Real> & aka_nodes = this->aka_mesh->getNodes();
Uint nb_nodes = this->nodes.getSize();
for (Uint node_id = 0; node_id < nb_nodes ; ++node_id) {
aka_nodes.push_back(&this->nodes.getVector().at(DIM*node_id));
}
// fill in the connectivity list of the aka_mesh
akantu::Array<Uint> * aka_connectivity;
if (DIM == twoD) {
this->aka_mesh->addConnectivityType(akantu::_triangle_3);
aka_connectivity = &this->aka_mesh->getConnectivity(akantu::_triangle_3);
} else if (DIM == threeD) {
this->aka_mesh->addConnectivityType(akantu::_tetrahedron_4);
aka_connectivity = &this->aka_mesh->getConnectivity(akantu::_tetrahedron_4);
}
Uint nb_elems = this->connectivity.size()/(DIM+1);
Uint elem[DIM+1];
for (Uint elem_id = 0; elem_id < nb_elems ; ++elem_id) {
// yeah, yeah, double copy, can't be bothered
for (Uint i = 0; i < DIM+1 ; ++i) {
elem[i] = static_cast<Uint>(this->connectivity.at((DIM+1)*elem_id + i));
}
if (signOfJacobian<DIM>(elem, this->nodes.getVector()) < 0) {
Uint tmp = elem[0];
elem[0] = elem[1];
elem[1] = tmp;
}
aka_connectivity->push_back(elem);
}
this->aka_mesh_exists = true;
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Mesh<DIM>::dump_msh(std::string filename) throw (std::string) {
if (!this->aka_mesh_exists) {
this->create_aka_mesh();
}
akantu::MeshIOMSH mesh_io;
mesh_io.write(filename, *this->aka_mesh);
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Mesh<DIM>::tagBoundaryPlane(Uint dir, Real val) {
if (!this->aka_mesh_exists) {
this->create_aka_mesh();
}
akantu::MeshUtils::buildFacets(*this->aka_mesh);
akantu::Array<Real> & nodes = this->aka_mesh->getNodes();
// find the elements of one dimension lower than the problem, they're the
// facets
typedef akantu::Array<Uint> conn_type;
akantu::Mesh::type_iterator surface_type = this->aka_mesh->firstType(DIM-1);
akantu::Mesh::type_iterator end_surface = this->aka_mesh->lastType(DIM-1);
for (; surface_type != end_surface; ++surface_type) {
conn_type & conn = this->aka_mesh->getConnectivity(*surface_type);
Uint nb_per_elem = this->aka_mesh->getNbNodesPerElement(*surface_type);
conn_type new_conn(0, nb_per_elem);
auto tag_array = this->aka_mesh->template getDataPointer<Uint>("tag_0", *surface_type);
tag_array->resize(0);
auto element = conn.begin(nb_per_elem);
auto end_element = conn.end(nb_per_elem);
for (; element != end_element; ++element) {
bool is_in_plane = true;
for (Uint node_id = 0 ; node_id < nb_per_elem ; ++node_id) {
is_in_plane &= (&nodes((*element)(node_id)))[dir] == val;
}
if (is_in_plane) {
tag_array->push_back(20);
new_conn.push_back(*element);
}
}
conn.resize(new_conn.getSize());
conn.copy(new_conn);
}
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Mesh<DIM>::tagBoundaryGeometry(const Geometry<DIM> & boundary_geom) {
if (!this->aka_mesh_exists) {
this->create_aka_mesh();
}
akantu::MeshUtils::buildFacets(*this->aka_mesh);
akantu::Array<Real> & nodes = this->aka_mesh->getNodes();
// find the elements of one dimension lower than the problem, they're the
// facets
typedef akantu::Array<Uint> conn_type;
akantu::Mesh::type_iterator surface_type = this->aka_mesh->firstType(DIM-1);
akantu::Mesh::type_iterator end_surface = this->aka_mesh->lastType(DIM-1);
for (; surface_type != end_surface; ++surface_type) {
conn_type & conn = this->aka_mesh->getConnectivity(*surface_type);
Uint nb_per_elem = this->aka_mesh->getNbNodesPerElement(*surface_type);
conn_type new_conn(0, nb_per_elem);
auto tag_array = this->aka_mesh->template getDataPointer<Uint>("tag_0", *surface_type);
tag_array->resize(0);
auto element = conn.begin(nb_per_elem);
auto end_element = conn.end(nb_per_elem);
for (; element != end_element; ++element) {
bool is_in_geom = true;
for (Uint node_id = 0 ; node_id < nb_per_elem ; ++node_id) {
is_in_geom &= boundary_geom.is_inside(&nodes((*element)(node_id)), 0);
}
if (is_in_geom) {
tag_array->push_back(20);
new_conn.push_back(*element);
}
}
conn.resize(new_conn.getSize());
conn.copy(new_conn);
}
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
inline bool compare_points(const PointRef<DIM> & point1,
const PointRef<DIM> & point2) {
bool retval = true;
for (Uint i = 0; i < DIM ; ++i) {
retval &= (fabs(1 - point1.getComponent(i)/point2.getComponent(i)) <
std::numeric_limits<Real>::epsilon());
}
return retval;
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Mesh<DIM>::compute_interface_atoms(PointContainer<DIM> & container,
const Geometry<DIM> & atomistic,
Real tol) {
if (!this->aka_mesh_exists) {
this->create_aka_mesh();
}
akantu::MeshUtils::buildFacets(*this->aka_mesh);
akantu::Array<Real> & nodes = this->aka_mesh->getNodes();
// get the set of all surface node ids
std::set<Uint> interface_point_ids;
this->surface_nodes(interface_point_ids);
typedef std::set<Uint>::iterator interface_iterator;
interface_iterator point_id = interface_point_ids.begin();
interface_iterator end_point_id = interface_point_ids.end();
for (; point_id != end_point_id ; ++point_id) {
const PointRef<DIM> point(&nodes(*point_id));
if (atomistic.is_inside(point, tol)) {
container.addPoint(point);
}
}
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Mesh<DIM>::surface_nodes(std::set<Uint> & indices) {
indices.clear();
if (!this->aka_mesh_exists) {
this->create_aka_mesh();
}
akantu::MeshUtils::buildFacets(*this->aka_mesh);
// find the elements of one dimension lower than the problem, they're the
// facets
typedef akantu::Array<Uint> conn_type;
akantu::Mesh::type_iterator surface_type = this->aka_mesh->firstType(DIM-1);
akantu::Mesh::type_iterator end_surface_type = this->aka_mesh->lastType(DIM-1);
for (; surface_type != end_surface_type; ++surface_type) {
const conn_type & conn = this->aka_mesh->getConnectivity(*surface_type);
auto element = conn.begin(DIM);
auto end_element = conn.end(DIM);
for (; element != end_element; ++element) {
for (Uint node = 0; node < DIM ; ++node) {
indices.insert((*element)(node));
}
}
}
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
std::string Mesh<DIM>::aka_mesh_name = "mesh 1";
template <Uint DIM>
Uint Mesh<DIM>::aka_mesh_counter = 1;
double invert(double x){return 1./x;}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
std::vector<Real> & Mesh<DIM>::computeRadiusRatios(Uint index)
{
Uint nb_elements;
int * tri_node;
if (index == std::numeric_limits<Uint>::max()) {
nb_elements = this->getNbElems();
tri_node = &this->connectivity.front();
this->radius_ratio_computed = true;
} else {
nb_elements = 1;
tri_node = &this->connectivity.at((DIM+1)*index);
this->radius_ratio_computed = false;
}
Real q_min= std::numeric_limits<Real>::max(), q_max=0;
this->radius_ratios.resize(nb_elements);
int nb_points = this->nodes.getSize();
const Real * points_z = &this->nodes.getVector().front();
int tri_order = DIM + 1; // it's weird, but this seems to be what it wants
int tri_num = static_cast<int>(nb_elements);
Real q_ave, q_area, q_var; //garbage, throwaway
Real * q_vec = &this->radius_ratios.front();
int offset = 0;
if (DIM == twoD) {
triangle_quality::q_measure(nb_points, points_z, tri_order, tri_num,
tri_node, &q_min, &q_max, &q_ave, &q_area,
q_vec, offset);
} else if (DIM == threeD) {
tetrahedron_quality::tet_mesh_quality1(nb_points, points_z, tri_order, tri_num,
tri_node, &q_min, &q_ave, &q_max, &q_var,
q_vec, offset);
} else {
std::stringstream error_stream;
error_stream << "computeRadiusRatios not implemented for DIM = " << DIM;
throw (error_stream.str());
}
std::vector<Real>::iterator ratio = this->radius_ratios.begin();
std::vector<Real>::iterator end_ratio = this->radius_ratios.end();
std::transform(ratio, end_ratio, ratio, invert);
this->radius_ratio_max = 1/q_min;
this->radius_ratio_min = 1/q_max;
return this->radius_ratios;
}
/* -------------------------------------------------------------------------- */
template<typename numtype>
class greater_than_object {
public:
greater_than_object(const numtype & threshold_):threshold(threshold_) {}
bool operator() (numtype i) {
return i > this->threshold;
}
private:
numtype threshold;
};
/* -------------------------------------------------------------------------- */
template<Uint DIM>
void Mesh<DIM>::getDegenerateElements(const Real & quality_threshold,
std::vector<int> & indices) {
if (!this->radius_ratio_computed) {
this->computeRadiusRatios();
}
greater_than_object<Real> comparator(quality_threshold);
std::vector<Real>::const_iterator begin_ratio = this->radius_ratios.begin();
std::vector<Real>::iterator ratio = this->radius_ratios.begin();
std::vector<Real>::iterator end_ratio = this->radius_ratios.end();
ratio = find_if(ratio, end_ratio, comparator);
while (ratio != end_ratio) {
indices.push_back(ratio - begin_ratio);
ratio = find_if(++ratio, end_ratio, comparator);
}
}
/* -------------------------------------------------------------------------- */
template<Uint DIM>
void Mesh<DIM>::getElementPlane(const int & elem_index, Real plane[DIM+1]) {
Real rhs[DIM+1], lhs[DIM*(DIM+1)];
int * element = &this->connectivity.at((DIM+1)*elem_index);
for (Uint node_id = 0 ; node_id < DIM+1 ; ++node_id) {
const Real * node_coord = &this->nodes.getVector().at(DIM*element[node_id]);
rhs[node_id] = node_coord[DIM-1];
Uint i;
for ( i = 0 ; i < DIM-1 ; ++i) {
lhs[DIM*node_id + i] = node_coord[i];
}
lhs[DIM*node_id + i] = 1;
}
// multiplication X'*X
Real XX[DIM*DIM];
cblas_dgemm(CblasRowMajor,// const enum CBLAS_ORDER Order,
CblasTrans, // const enum CBLAS_TRANSPOSE TransA,
CblasNoTrans, // const enum CBLAS_TRANSPOSE TransB,
DIM, // const int M,
DIM, // const int N,
DIM+1, // const int K,
1., // const double alpha,
lhs, // const double *A,
DIM, // const int lda,
lhs, // const double *B,
DIM, // const int ldb,
0., // const double beta,
XX, // double *C,
DIM); // const int ldc);
//invert XX
inverse(XX, DIM);
// compute inv(X'X)^*X'
Real XXX[DIM*(DIM+1)];
cblas_dgemm(CblasRowMajor, // const enum CBLAS_ORDER Order,
CblasNoTrans, // const enum CBLAS_TRANSPOSE TransA,
CblasTrans, // const enum CBLAS_TRANSPOSE TransB,
DIM, // const int M,1
DIM+1, // const int N,
DIM, // const int K,
1., // const double alpha,
XX, // double *A,
DIM, // const int lda,
lhs, // const double *B,
DIM, // const int ldb,
0., // const double beta,
XXX, // double *C,
DIM + 1); // const int ldc);
/// compute tilde a
Real tilde_a[DIM] = {0};
cblas_dgemv(CblasRowMajor, // const enum CBLAS_ORDER Order,
CblasNoTrans, // const enum CBLAS_TRANSPOSE TransA,
DIM, // const int M,
DIM + 1, // const int N,
1., // const double alpha,
XXX, // const double *A,
DIM + 1, // const int lda,
rhs, // const double *X,
1, // const int incX,
0., // const double beta,
tilde_a, // double *Y,
1); // const int incY);
Real denumerator = 1;
for (Uint j = 0 ; j < DIM-1 ; ++j) {
denumerator += tilde_a[j]*tilde_a[j];
}
Real c = pow(1/denumerator, 0.5);
Uint j;
for ( j = 0 ; j < DIM-1 ; ++j) {
plane[j] = tilde_a[j]*c;
}
plane[j] = c;
plane[j+1] = tilde_a[j]*c;
}
/* -------------------------------------------------------------------------- */
template<Uint DIM>
void Mesh<DIM>::getCircumSphere(const int & elem_index, Real& radius, Real centre [DIM]) {
Real tetra[DIM*(DIM+1)];
int * element = &this->connectivity.at((DIM+1)*elem_index);
for (Uint node = 0 ; node < DIM+1 ; ++node) {
const Real * node_coord = &this->nodes.getVector().at(DIM*element[node]);
for (Uint dir = 0 ; dir < DIM ; ++dir) {
tetra[node*DIM + dir] = node_coord[dir];
}
}
if (DIM == twoD) {
// compute the lengths of the triangle
Real lengths[DIM+1];
for (Uint vert_id = 0 ; vert_id < DIM+1 ; ++vert_id) {
lengths[vert_id] = 0;
for (Uint dir = 0 ; dir < DIM ; ++dir) {
lengths[vert_id] += square(tetra[DIM*vert_id+dir] -
tetra[DIM*((vert_id+1)%(DIM+1))+dir]);
}
lengths[vert_id] = sqrt(lengths[vert_id]);
}
Real semi_perim = 0;
for (Uint i = 0 ; i < DIM+1 ; ++i) {
semi_perim += 0.5*lengths[i];
}
radius = lengths[0]*lengths[1]*lengths[2]/
(4*sqrt(semi_perim *
(semi_perim - lengths[0]) *
(semi_perim - lengths[1]) *
(semi_perim - lengths[2])));
Real D_param = 0;
for (Uint vert_id = 0 ; vert_id < 1+DIM ; ++vert_id) {
D_param += 2*(tetra[DIM*vert_id] *
(tetra[DIM*((vert_id+1)%(DIM+1))+1] -
tetra[DIM*((vert_id+2)%(DIM+1))+1]));
}
for (Uint dir = 0 ; dir < DIM ; ++dir) {
centre[dir] = 0;
for (Uint vert_id = 0 ; vert_id < DIM+1 ; ++vert_id) {
Uint a_id = DIM*vert_id;
Uint b_id = DIM*((vert_id+1)%(DIM+1)) + (dir+1)%2;
Uint c_id = DIM*((vert_id+2)%(DIM+1)) + (dir+1)%2;
centre[dir] += dot_prod<DIM>(&tetra[a_id], &tetra[a_id]) *
(tetra[b_id]-tetra[c_id])/D_param;
}
}
} else if (DIM == threeD) {
tetrahedron_quality::tetrahedron_circumsphere_3d(tetra, &radius, centre);
} else {
throw (std::string("getCircumSphere has not been implemented for this dimension"));
}
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Mesh<DIM>::splitDegenerates(const Real & radius_threshold,
PointContainer<DIM> & points) {
typedef std::vector<int> intvec;
intvec indices;
this->getDegenerateElements(radius_threshold, indices);
intvec::iterator elem = indices.begin();
intvec::iterator end_elem = indices.end();
for (; elem != end_elem; ++elem) {
Real point[DIM] = {0};
const int * node_ids = this->getElem(*elem);//&this->connectivity.at((DIM+1)*(*elem));
for (Uint node = 0 ; node < DIM+1 ; ++node) {
const Real * coord = &this->nodes.getVector().at(DIM*node_ids[node]);
for (Uint dir = 0 ; dir < DIM ; ++dir) {
point[dir] += coord[dir]/(DIM+1);
}
}
points.addPoint(point);
}
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Mesh<DIM>::get_neighbour_lists(std::vector<std::set<int>> & neighbours) {
Uint nb_nodes = this->nodes.getSize();
neighbours.clear();
neighbours.resize(nb_nodes);
Uint nb_elems = this->getNbElems();
for (Uint elem_id = 0 ; elem_id < nb_elems ; ++elem_id) {
const int * elem = this->getElem(elem_id);
for (Uint local_nd_nb = 0 ; local_nd_nb < DIM+1 ; ++local_nd_nb) {
const Uint global_nd_id = elem[local_nd_nb];
neighbours.at(global_nd_id).insert(elem_id);
}
}
Uint maxneigh = 0, max_id = std::numeric_limits<Uint>::max();
for (Uint i = 0 ; i < neighbours.size() ; ++i) {
Uint nb_neigh = neighbours[i].size();
if (nb_neigh > maxneigh) {
maxneigh = nb_neigh;
max_id = i;
}
}
std::cout << "Element " << max_id << " has the maximum number of neighbours with "
<< maxneigh << " neighbours." <<std::endl;
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Mesh<DIM>::extract_neighbours(Uint elem_id,
const std::vector<std::set<int> > & neighbours,
const std::set<Uint> surf_ids,
const Geometry<DIM> & atomic,
std::set<int> & extracted_elem_ids,
std::set<int> & free_nodes,
std::set<Uint> & degen_ids) {
extracted_elem_ids.insert(elem_id);
degen_ids.erase(degen_ids.find(elem_id));
const int * elem = this->getElem(elem_id);
for (Uint i = 0 ; i < DIM+1 ; ++i) {
bool is_not_surface_nd = (surf_ids.find(elem[i]) == surf_ids.end());
PointRef<DIM> point(this->nodes[elem[i]]);
bool is_not_atom = !atomic.is_inside(point, 0);
if (is_not_surface_nd && is_not_atom) {
free_nodes.insert(elem[i]);
}
for (Uint neigh_id: neighbours.at(elem[i])) {
bool is_new = (extracted_elem_ids.find(neigh_id) == extracted_elem_ids.end());
if (is_new) {
extracted_elem_ids.insert(neigh_id);
bool is_degenerate = (degen_ids.find(neigh_id) != degen_ids.end());
if (is_degenerate) {
this->extract_neighbours(neigh_id, neighbours, surf_ids, atomic,
extracted_elem_ids, free_nodes, degen_ids);
}
}
}
}
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Mesh<DIM>::relaxDegenerates(const Real & radius_threshold,
const Geometry<DIM> & atomic,
const Real & step_size,
const Real & min_step_size,
const Uint & max_iter,
bool blocked_nodes_lethal,
bool verbose) {
std::set<Uint> surf_nd_ids, degen_ids;
// get the set of surface nodes (they are always blocked)
this->surface_nodes(surf_nd_ids);
typedef std::vector<int> intvec;
intvec indices;
this->getDegenerateElements(radius_threshold, indices);
// get a set of degenerate elements (for easier erasing)
for (const auto & id: indices) {
degen_ids.insert(id);
}
std::vector<std::set<int>> neighbours;
this->get_neighbour_lists(neighbours);
// loop through all the degenerates
while (!degen_ids.empty()) {
if (verbose) {
std::cout<< "There are " << degen_ids.size()
<< " degenerates left to deal with" << std::endl;
}
std::set<int> extracted_elem_ids;
std::set<int> degenerate_neighbours;
std::set<int> free_nodes;
int current_elem_id = *degen_ids.begin();
if (verbose) {
std::cout << "relaxing element " << current_elem_id
<< " with quality threshold " << radius_threshold
<< ". current quality is "
<< computeRadiusRatio<DIM>(this->nodes.getArray(),
this->getElem(current_elem_id),
this->nodes.getSize())
<< std::endl;
}
this->extract_neighbours(current_elem_id,
neighbours,
surf_nd_ids,
atomic,
extracted_elem_ids,
free_nodes,
degen_ids);
if (free_nodes.empty()) {
std::stringstream error, name;
name << "nodes of element " << current_elem_id;
PointContainer<DIM> elem_nodes(name.str());
const int * elem = this->getElem(current_elem_id);
for (Uint i = 0 ; i < DIM+1 ; ++i) {
elem_nodes.addPoint(this->nodes[elem[i]]);
}
Real quality = computeRadiusRatio<DIM>(this->nodes.getArray(),
elem, this->nodes.getSize());
error << "There are no free nodes in element " << current_elem_id
<< ". The nodes are " << std::endl << elem_nodes << std::endl
<< ", and the quality is " << quality
<< ". Make sure the element is not in the pad or spanning through "
<< "the full thickness of the domain." <<std::endl;
if (blocked_nodes_lethal) {
throw std::runtime_error(error.str());
} else {
std::cout << "WARNING! " << error.str();
continue;
}
}
//for each extracted neighbourhood, relax the degenerates
std::vector<const int*> neighbourhood_connectivity;
for (const int & elem_id: extracted_elem_ids) {
neighbourhood_connectivity.push_back(this->getElem(elem_id));
}
std::vector<int> free_nodes_vec;
for (const int & node_id: free_nodes) {
free_nodes_vec.push_back(node_id);
}
MiniMesh<DIM> neighbourhood(neighbourhood_connectivity,
free_nodes_vec,
this->nodes);
try {
neighbourhood.relax(step_size, max_iter, radius_threshold, min_step_size,
verbose);
} catch (typename MiniMesh<DIM>::StopIterError & e) {
if (blocked_nodes_lethal) {
throw e;
} else {
std::cout << "WARNING! Element " <<current_elem_id << ": "
<< e.what() << std::endl;
}
}
}
}
/* -------------------------------------------------------------------------- */
template class Mesh<2>;
template class Mesh<3>;
c4science · Help