Homec4science
Diffusion CADD_mesher (dev/incorporate-qhull) Browse

rCADDMESH/src/domain/mesh.cc
6c4fc7e2918adev/incorporate-qhull

Last Change
Raw File

mesh.cc
View Options

extern "C" {
#include <stdio.h>
#include <qhull/qhull.h>
#include <qhull/mem.h>
#include <qhull/qset.h>
#include <qhull/poly.h> /* for qh_vertexneighbors() */
#include <qhull/geom.h> /* for qh_facetarea(facet)*/
}
#include "mesh.hh"
#include <cblas.h>
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);
// general matrix vector multiplication
}
template <Uint DIM>
Mesh<DIM>::Mesh(const std::vector<Real> & points, Geometry<DIM> * excluded_points)
{
for (Uint i = 0; i < points.size()/DIM ; i += DIM) {
if (excluded_points.is_inside(&points.at(i)) <= 0 ) {
for (Uint j = 0; j < DIM ; ++j) {
this->nodes.push_back(points.at(i+j));
}
}
}
}
template <Uint DIM>
Mesh<DIM>::printself(std::ostream & stream, int indent) const
{
indent_space(stream, indent);
stream << "PRINTSELF FOR MESH TO BE IMPLEMENTED" << std::endl;
}
template <Uint DIM>
inline bool degenerate_finder(const std::vector<double> & points,
const std::vector<int> & element,
double volume,
double criterion)
{
double dist = 0;
for (unsigned int i = 0; i < DIM; ++i) {
double ny_dist = fabs(points[DIM*element[0]+i]-points[DIM*element[1]+i]);
keep_max(dist, ny_dist);
}
std::cout << pow(dist, DIM)/(DIM*criterion) << "\t>\t" << volume << "\t";
return volume > pow(dist, DIM)/(DIM*criterion);
}
template <Uint DIM>
Mesh<DIM>::delaunay(Real degeneration_criterion){
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 */
/* Call qhull */
exitcode = qh_new_qhull(dim, this->nodes.size()/DIM, &this->nodes[0],
ismalloc, flags, NULL, stderr);
if (exitcode){
std::stringstream exit_stream;
exit_stream << "qhull failed with exitcode " << exitcode;
throw exit_stream.str();
} else {
//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 << "Element #" << nb_elems
<< "has wrong number of vertices: " << counter;
throw exit_stream.str();
} else {
bool is_valid = degenerate_finder(this->nodes, current_element,
DIM, element_volume, degeneration_criterion);
if (is_valid) {
for (unsigned int i = 0; i < DIM+1; ++i) {
connectivity.push_back(current_element[i]);
}
}
}
}
}
/* Free the memory */
qh_freeqhull(qh_ALL);
}
template <Uint DIM>
Mesh<DIM>cleanup(const Geometry & atomistic_domain) {
Uint nb_points = this->nodes.size()/DIM;
Uint nb_elems = this->connectivity.size()/DIM;
std::vector<Real> new_points;
std::vector<int> new_connectivity;
std::vector<int> permutation(nb_points, -1);
int current_element[DIM+1];
Uint counter = 0;
// loop through elements
for (Uint i = 0; i < nb_elems ; ++i) {
int * element_nodes &connectivity.at(i*(DIM+1));
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];
if (atomistic_domain.is_inside(&nodes.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) {
for (Uint j = 0; j < DIM+1 ; ++j) {
int point_id = element_nodes[j];
if (permutation.at(point_id) < 0) {
permutation.at(point_id) = counter;
++counter;
for (Uint k = 0; k < DIM ; ++k) {
new_points.push_back(this->nodes.at(point_id*DIM+k));
}
}
new_connectivity.push_back(permutation.at(point_id));
}
}
}
this->nodes = new_points;
this->connectivity = new_connectivity;
}
void inverse(double * A, int N)
{
int *IPIV = new int[N+1];
int LWORK = N*N;
double *WORK = new double[LWORK];
int INFO;
dgetrf_(&N,&N,A,&N,IPIV,&INFO);
dgetri_(&N,A,&N,IPIV,WORK,&LWORK,&INFO);
delete IPIV;
delete WORK;
}
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)));
Real transformation[(DIM+1)*(DIM+1)] = {0};
// 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.at(element_nodes[DIM]*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.at(element_nodes[node_local]*DIM + dir)
- last_node_coord);
}
}
//invert the matrix to make it useful
inverse(&(this->reverse_matrices.at(elem_id*(DIM+1)*(DIM+1))), DIM + 1);
}
}
template<Uint DIM>
Real Mesh<DIM>::interpolate(Real * coords) const {
Uint nb_elems = this->connectivity.size()/(DIM+1);
Real elem_coords[DIM+2]; //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];
}
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) {
Real* matrix = &(this->reverse_matrices.at(elem_id * (DIM+1)*(DIM+1)));
cblas_dgemv(CblasRowMajor, CblasNoTrans, DIM+1, DIM+1, 1., 0., matrix,
glob_coords, elem_coords, DIM+1, 1, 1);
if (fabs(1-elem_coords[0]) > std::numeric_limits<Real>::epsilon()){
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]);
if (minimum >= 0) {
found = true;
element_nodes = &(this->connectivity.at(i*(DIM+1)));
break;
}
}
if (!found) {
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] ") are not within the domain. " << std::endl;
throw(error_stream.str());
}
Real interpolated_value = 0;
for (Uint i = 0; i < DIM+1 ; ++i) {
interpolated_value += this->nodal_vals.at(element_nodes[i]) * shape_vals[i];
}
return interpolated_value;
}
c4science · Help