diff --git a/src/domain/cadd_mesh.cc b/src/domain/cadd_mesh.cc
index e2c7376..6bfa912 100644
--- a/src/domain/cadd_mesh.cc
+++ b/src/domain/cadd_mesh.cc
@@ -1,796 +1,830 @@
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 "cadd_mesh.hh"
#include <cmath>
#include <sstream>
#include <limits>
#include <set>
#include <algorithm>
#include <stdexcept>
#include <cblas.h>
#include "dumper_paraview.hh"
// akantu includes for mesh io
#include <aka_common.hh>
#include <aka_vector.hh>
#include <mesh_io_msh.hh>
#include <mesh_utils.hh>
// external routines for quality measures
#include "tet_mesh_quality.hh"
#include "triangulation_quality.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, 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)
{
for (Uint i = 0; i < points.size() ; i += DIM) {
if ((excluded_points == NULL) || (!excluded_points->is_inside(&points.at(i))) ) {
this->renumerotation.push_back(i/DIM);
for (Uint j = 0; j < DIM ; ++j) {
this->nodes.push_back(points.at(i+j));
}
}
}
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, "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 Real computeRadiusRatio(Real * points, int * element, Uint nb_points) {
Uint nb_elements = 1;
Real q_min= std::numeric_limits<Real>::max(), q_max=0;
int tri_order = DIM + 1; // it's weird, but this seems to be what it wants
Real q_ave, q_area, q_var; //garbage, throwaway
Real ratio = 0;
int offset = 0;
if (DIM == twoD) {
triangle_quality::q_measure(nb_points, points, tri_order, nb_elements,
element, &q_min, &q_max, &q_ave, &q_area,
&ratio, offset);
} else if (DIM == threeD) {
tetrahedron_quality::tet_mesh_quality1(nb_points, points, tri_order, nb_elements,
element, &q_min, &q_ave, &q_max, &q_var,
&ratio, offset);
} else {
std::stringstream error_stream;
error_stream << "computeRadiusRatio not implemented for DIM = " << DIM;
throw (error_stream.str());
}
return 1./ratio;
}
template <Uint DIM>
inline bool degenerate_finder(std::vector<double> & points,
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);
//}
//return volume > pow(dist, DIM)/(DIM*criterion);
Uint nb_nodes = points.size()/DIM;
Real * points_ptr = &points.front();
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>
void Mesh<DIM>::delaunay(const Real & degeneration_criterion){
// make a copy of the nodes before they get jiggled
std::vector<Real> copy_nodes = this->nodes;
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, copy_nodes.size()/DIM, ©_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<DIM>(this->nodes, current_element,
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);
this->meshed = true;
}
template <Uint DIM>
void Mesh<DIM>::cleanup(const Geometry<DIM> & atomistic_domain) {
Uint nb_points = this->nodes.size()/DIM;
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] += nodes.at(DIM*point_id + k);
}
if (atomistic_domain.is_inside(&nodes.at(DIM*point_id)) ) {
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)) {
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.at(point_id*DIM+k));
}
}
new_connectivity.push_back(this->permutation->at(point_id));
}
}
}
}
this->nodes = 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.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);
}
this->interpolation_initialised = true;
}
template<Uint DIM>
Real Mesh<DIM>::interpolate(const Real * coords,
const std::vector<Real> & nodal_values) 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];
}
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]) > std::numeric_limits<Real>::epsilon()){
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 + 2*std::numeric_limits<Real>::epsilon() >= 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 = " << std::numeric_limits<Real>::epsilon()
<< 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];
}
index = this->renumerotation.at(index);
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);
dumper.setPoints(&this->nodes[0], DIM, this->nodes.size()/DIM, 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.size()/DIM;
for (Uint node_id = 0; node_id < nb_nodes ; ++node_id) {
aka_nodes.push_back(&this->nodes.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) < 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) {
+ const conn_type & conn = this->aka_mesh->getConnectivity(*surface_type);
+
+ auto tag_array = this->aka_mesh->getUIntDataPointer(*surface_type, "tag_1");
+
+ auto element = conn.begin(DIM);
+ auto end_element = conn.end(DIM);
+ for (; element != end_element; ++element) {
+ bool is_in_plane = true;
+ for (Uint node_id = 0 ; node_id < DIM ; ++node_id) {
+ is_in_plane &= (&nodes((*element)(node_id)))[dir] == val;
+ }
+ if (is_in_plane) {
+ tag_array->push_back(2);
+ } else {
+ tag_array->push_back(1);
+ }
+ }
+ }
+}
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) {
if (!this->aka_mesh_exists) {
this->create_aka_mesh();
}
akantu::MeshUtils::buildFacets(*this->aka_mesh);
akantu::Array<Real> & nodes = this->aka_mesh->getNodes();
std::set<Uint> interface_point_ids;
// 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) {
const conn_type & conn = this->aka_mesh->getConnectivity(*surface_type);
typedef akantu::Array<Uint>::const_iterator
<akantu::Vector<Uint> > element_iterator;
element_iterator element = conn.begin(DIM);
element_iterator end_element = conn.end(DIM);
for (; element != end_element; ++element) {
for (Uint node = 0; node < DIM ; ++node) {
PointRef<DIM> point(&nodes((*element)(node)));
if (atomistic.is_inside(point)) {
interface_point_ids.insert((*element)(node));
}
}
}
}
//the set is now filled with the indices of all interface atoms
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));
container.addPoint(point);
}
}
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->connectivity.size()/(DIM+1);
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.size()/DIM;
Real * points_z = &this->nodes.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) {
Real * node_coord = &this->nodes.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) {
Real * node_coord = &this->nodes.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};
int * node_ids = &this->connectivity.at((DIM+1)*(*elem));
for (Uint node = 0 ; node < DIM+1 ; ++node) {
Real * coord = &this->nodes.at(DIM*node_ids[node]);
for (Uint dir = 0 ; dir < DIM ; ++dir) {
point[dir] += coord[dir]/(DIM+1);
}
}
points.addPoint(point);
}
}
template class Mesh<2>;
template class Mesh<3>;
diff --git a/src/domain/cadd_mesh.hh b/src/domain/cadd_mesh.hh
index 7932abd..3d16c8d 100644
--- a/src/domain/cadd_mesh.hh
+++ b/src/domain/cadd_mesh.hh
@@ -1,219 +1,228 @@
/**
* @file mesh.hh
* @author Till Junge <junge@lsmspc42.epfl.ch>
* @date Fri Apr 20 14:52:26 2012
*
* @brief
*
* @section LICENSE
*
* <insert lisence here>
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __CADD_MESHER_MESH_HH__
#define __CADD_MESHER_MESH_HH__
#include "common.hh"
#include "geometry.hh"
#include "mesh.hh" //this is the akantu mesh class
#include <vector>
/*! \class Mesh builds FEM mesh
*/
template <Uint DIM>
class Mesh {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
/*! Constructor
* \param points vector (e.g., from a PointContainer) with the nodal points to
* mesh
* \param excluded_points Geometry which contains points to be excluded from
* the mesh (e.g., the pure atoms (non-pad, non-interface) in the case of
* CADD). This also allows for non-convex geometries.
*/
Mesh(const std::vector<Real> & points, Geometry<DIM> * excluded_points = NULL);
virtual ~Mesh();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// function to print the contain of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
/*! meshes it all up. the degeneration_criterion d is used to evaluate whether
* a element is a degenerate sliver. 1000 seems to be an ok threshold for this
* \param degeneration_criterion
*/
void delaunay(const Real & degeneration_criterion);
/*! excludes all elements whose nodes are \a all within \a atomistic_domain
* or whose centre of gravity is whithin \a atomistic_domain
* \param atomistic_domain Geometry instance*/
void cleanup(const Geometry<DIM> & atomistic_domain);
/*! for each element, compute the transformation matrix to obtain nodal
* coordinates from global carthesian ones:
* for 2D the shape functions are:
*\f[N_1 = \xi, N_2 = \eta, N_3 = 1-\xi-\eta\mbox{\ \ \ }(1)\f]
* They can be used to interpolate nodal coordinates \f$\boldsymbol\xi = (\xi, \eta)\f$ to obtain global coordinates \f$\boldsymbol x = (x, y)\f$:
* \f[\boldsymbol x = \sum_{i=1}^3 N_1(\boldsymbol\xi)\boldsymbol x_i\mbox{\ \ \ }(2)\f]
Substituting (1) into (2) yields the system of equations
\f[
\left(\begin{array}{c}x\\y\end{array}\right) = \left(\begin{array}{cc}
x_1-x_3 & x_2-x_3 \\
y_1-y_3 & y_2-y_3 \\
\end{array}\right)\left(\begin{array}{c}\xi\\\eta\end{array}\right)
+ \left(\begin{array}{c}x_3\\y_3\end{array}\right)
\f]
I've homogenised the coordinates to incorporate the translation part:
\f[
\boldsymbol x = \left(\begin{array}{c}1\\x\\y\end{array}\right)=
\underbrace{\left(\begin{array}{ccc}
1 & 0 & 0 \\
x_3 & x_1-x_3 & x_2-x_3 \\
y_3 & y_1-y_3 & y_2-y_3 \\
\end{array}\right)}_{\mathbf{S}}
\left(\begin{array}{c}1\\\xi\\\eta\end{array}\right)
\f]
This function computes the inverse of \f$\mathbf S\f$ for each element
*/
void initInterpolation();
/*!interpolate \a nodal_values at the point \a coords. This method can only be
* used after initInterpolation() has been called
* \param coords
* \param nodal_values These nodal values have to be in the order of the
* original \a points vector given to the constructor.
*/
Real interpolate(const Real * coords, const std::vector<Real> & nodal_values) const;
inline Real interpolate(const PointRef<DIM> & point,
const std::vector<Real> & nodal_values) const{
return this->interpolate(point.getArray(), nodal_values);
}
/*!writes a paraview file containing the nodes and
* connectivities for visualisation
\param filename path of the output file, starting from "./"*/
void dump_paraview(std::string filename) throw (std::string) ;
/*!writes a msh file containing the nodes and
* connectivities for visualisation
\param filename path of the output file, starting from "./"*/
void dump_msh(std::string filename) throw (std::string) ;
+ /*! Compute the surface elements of one boundary surface and tag them
+ * needed to allow libmultiscale to apply natural boundary conditions
+ * surface is given by the direction and offset of the plane it's in.
+ * only planes normal to a cooridnate axis possible for now
+ * the boundary surface is tagged as number 2
+ \param dir dimension in which the normal to the plane points
+ \param val offset
+ */
+ void tagBoundaryPlane(Uint dir, Real val);
/*! Fills a PointContainer with interface atoms
\param container this will be filled with inderface atoms
\param atomistic domain in which interface atoms are expected
*/
void compute_interface_atoms(PointContainer<DIM> & container,
const Geometry<DIM> & atomistic);
/*! conputes the radius ratio for every element
*/
std::vector<Real> & computeRadiusRatios(Uint index = std::numeric_limits<Uint>::max());
/*! returns a vector of element \a indices for which the radius radio exceeds
* \a quality_threshold
* \param quality_threshold a positive real number.
* \param indices empty vector of indices
*/
void getDegenerateElements(const Real & quality_threshold,
std::vector<int> & indices);
/*! computes the plane (3d)/line (2d) which best approximates the positions
* of the elemnt's nodes:
* \f[ ax + by + cz + d = 0\f]
* \f[ \Rightarrow z = \underbrace{-\frac{a}{c}}_{\tilde a}x
* \underbrace{-\frac{b}{c}}_{\tilde b}y \underbrace{-\frac{d}{c}}_{\tilde d}\f]
* by defining \f$\boldsymbol z = \left(z_1, z_2, z_3, z_4\right)^\mathrm{T}\f$,
* \f$\mathbf{X} = \left[\begin{array}{ccc} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ x_3 & y_3 & 1 \\ x_4 & y_4 & 1\end{array}\right]\f$ and \f$\boldsymbol{\tilde a} = \left(\tilde a, \tilde b, \tilde d\right)^\mathrm{T}\f$ , we get
* \f[\boldsymbol{\tilde a} = \left(\mathbf{X} ^\mathrm{T}\mathbf{X}\right)^\mathrm{-1}\mathbf{X}^\mathrm{T}\boldsymbol z.\f]
* we enforce that the normal vector \f$\left|\left|\left(a,b,c\right)^\mathrm{T}\right|\right| = 1\f$ by setting by computing \f$c\f$ through \f[ \left(\frac{a}{c}\right)^2 + \left(\frac{b}{c}\right)^2 + 1 = \frac{1}{c^2}\f]
*/
void getElementPlane(const int & elem_index, Real plane[DIM+1]);
/*! compute the circumsphere (3d)/circumcircle (2d) of elemnt \a elem_index
* \param elem_index index of element
* \param radius the radius of the circumsphere
* \param centre coords of the center
*/
void getCircumSphere(const int & elem_index, Real& radius, Real centre [DIM]);
void splitDegenerates(const Real & radius_threshold,
PointContainer<DIM> & points);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/*! returns the vector of nodes*/
const std::vector<Real> & getNodes() const {return this->nodes;}
/*! returns the vector of connectivities*/
const std::vector<int> & getConnectivity() const {return this->connectivity;}
const Real & getRadiusRatioMax() {
if (!this->radius_ratio_computed) {
this->computeRadiusRatios();
}
return this->radius_ratio_max;
}
const Real & getRadiusRatioMin() {
if (!this->radius_ratio_computed) {
this->computeRadiusRatios();
}
return this->radius_ratio_min;
}
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
void create_aka_mesh();
// contains all the nodes
std::vector<Real> nodes;
// element definitions
std::vector<int> connectivity;
// transformation matrices stored by matrix to evaluate interpolated nodal
// values
std::vector<Real> reverse_matrices;
std::vector<int> * permutation, * remutation;
std::vector<int> renumerotation;
std::vector<Real> radius_ratios;
Real radius_ratio_max, radius_ratio_min;
bool radius_ratio_computed;
//
bool meshed, cleaned, interpolation_initialised;
akantu::Mesh * aka_mesh;
static std::string aka_mesh_name;
std::string my_aka_mesh_name;
static Uint aka_mesh_counter;
bool aka_mesh_exists;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
/// standard output stream operator
template <Uint DIM>
inline std::ostream & operator <<(std::ostream & stream, const Mesh<DIM> & _this)
{
_this.printself(stream);
return stream;
}
#endif /* __CADD_MESHER_MESH_HH__ */
diff --git a/src/domain/domain.cc b/src/domain/domain.cc
index 9513803..02bb9ca 100644
--- a/src/domain/domain.cc
+++ b/src/domain/domain.cc
@@ -1,255 +1,255 @@
#include "domain.hh"
#include "tree.hh"
#include <limits>
#include <cmath>
#include <iterator>
template <Uint DIM>
Domain<DIM>::Domain(const Geometry<DIM> & overall_,
const Geometry<DIM> & refined_,
const Geometry<DIM> & atomistic_,
const Geometry<DIM> & atomistic_skinned_,
const Lattice<DIM> & lattice_,
bool boundaries_special_treatment_,
const PointRef<DIM> & centre)
:auxiliary_mesh(NULL), tree(NULL), main_mesh(NULL),
filled(false), main_meshed(false),
boundaries_special_treatment(boundaries_special_treatment_),
mesh_quality(-std::numeric_limits<Real>::max())
{
this->overall = overall_.resolveType();
this->refined = refined_.resolveType();
this->atomistic = atomistic_.resolveType();
this->atomistic_skinned = atomistic_skinned_.resolveType();
this->lattice = lattice_.resolveType();
Real max_size = -std::numeric_limits<Real>::max();
for (Uint dir = 0; dir < DIM ; ++dir) {
Real size = this->overall->size_in_direction
(this->lattice->getLatticeVector(dir),
this->lattice->getConstant(dir));
keep_max(max_size, size);
}
// compute the next larger power of 2
int exponent = log(max_size)/log(2) +1;
Uint pow2size = pow(2, exponent);
if (¢re == NULL) {
this->tree = new Tree<DIM>(NULL, 2*pow2size, this->lattice, this);
} else {
this->tree = new Tree<DIM>(centre.getArray(), 2*pow2size, this->lattice, this);
}
}
template <Uint DIM>
Domain<DIM>::~Domain() {
delete this->auxiliary_mesh;
delete this->main_mesh;
delete this->tree;
}
template <Uint DIM>
void Domain<DIM>::setPointRefinement(const Real * coords, Real & refinement ) {
this->refinement_points.addPoint(coords);
this->refinement.push_back(refinement);
}
template <Uint DIM>
void Domain<DIM>::fill_points(bool periodicity[DIM]) throw(std::string){
this->initialiseAuxiliaryMesh();
if (this->boundaries_special_treatment) {
this->overall->generate_surface_points(*this->auxiliary_mesh,
this->refinement,
this->lattice->getConstants(),
periodicity,
this->points);
}
this->tree->fill(this->points);
this->fill_atoms();
this->filled = true;
this->points.clean_duplicates(1);
}
template<Uint DIM>
void Domain<DIM>::complement_periodically(int direction, Uint dim) throw(std::string){
if (!this->filled) {
throw("Mesh has got to be filled before it can be complemented");
}
std::cout << "Shit I've been called" << std::endl;
this->overall->complement_periodically(*this->auxiliary_mesh,
this->refinement,
this->points,
direction,
dim);
}
template <Uint DIM>
void Domain<DIM>::printself(std::ostream & stream, int indent) const{
indent_space(stream, indent);
stream << "PRINTSELF IS NOT YET IMOPLEMENTED FOR CLASS DOMAIN";
}
template <Uint DIM>
void Domain<DIM>::initialiseAuxiliaryMesh() {
this->auxiliary_mesh = new Mesh<DIM>(this->refinement_points.getVector());
this->auxiliary_mesh->delaunay(1000);
this->auxiliary_mesh->initInterpolation();
}
template <Uint DIM>
Mesh<DIM> & Domain<DIM>::getMesh(Real quality) throw(std::string) {
if ((!this->main_meshed) || (quality!=this->mesh_quality)) {
this->main_mesh = new Mesh<DIM>(this->points.getVector(), this->atomistic_skinned);
try {
this->main_mesh->delaunay(quality);
} catch (std::string & error) {
this->tree->dump();
this->points.dump_paraview();
throw("affen_fehler " + error);
}
this->main_mesh->cleanup(*this->atomistic);
this->main_meshed = true;
}
return *this->main_mesh;
}
template <Uint DIM>
void Domain<DIM>::splitDegenerates(const Real & radius_threshold) {
this->main_mesh->splitDegenerates(radius_threshold, this->points);
delete this->main_mesh;
this->main_meshed = false;
}
template <Uint DIM>
Mesh<DIM> & Domain<DIM>::getAuxiliaryMesh() throw(std::string) {
if (this->auxiliary_mesh == NULL) {
throw (std::string("Auxiliary mesh does not exist"));
}
return *this->auxiliary_mesh;
}
template <Uint DIM>
Tree<DIM> & Domain<DIM>::getTree() {
return *this->tree;
}
template <Uint DIM>
bool Domain<DIM>::inFemDomain(PointRef<DIM> & point) {
Real distance_from_border = -this->overall->from_border(point);
if (distance_from_border < 0) {
// we're out of the domain
return false;
}
Real refinement = int(this->interpolate(point))
* this->lattice->getRepresentativeConstant();
bool retval;
if (this->boundaries_special_treatment) {
retval = distance_from_border >= refinement
&& !this->atomistic_skinned->is_inside(point);
} else {
retval = !this->atomistic_skinned->is_inside(point);
}
return retval;
}
template <Uint DIM>
void Domain<DIM>::fill_atoms() {
Uint current_index = 0;
for (Uint point_id = 0; point_id < this->points.getSize() ; ++point_id) {
PointRef<DIM> point = this->points.getPoint(point_id);
if (this->refined->is_inside(point)) {
this->atoms.addPoint(point);
if (this->inFemDomain(point)) {
this->pad_atoms_indices.push_front(current_index);
}
++ current_index;
}
}
}
template <Uint DIM>
void Domain<DIM>::dump_atoms_lammps(std::string filename) {
if (this->atoms.getSize() == 0) {
this->fill_atoms();
}
this->atoms.dump_lammps(filename);
}
template <Uint DIM>
void Domain<DIM>::dump_atoms_paraview(std::string filename) {
if (this->atoms.getSize() == 0) {
this->fill_atoms();
}
this->atoms.dump_paraview(filename);
}
template <Uint DIM>
const PointContainer<DIM> & Domain<DIM>::getAtoms() {
if (this->atoms.getSize() == 0) {
this->fill_atoms();
}
return this->atoms;
}
template <Uint DIM>
bool Domain<DIM>::inDomain(PointRef<DIM> & point) {
Real distance_from_border = -this->overall->from_border(point);
if (distance_from_border < 0) {
// we're out of the domain
return false;
} else {
if (this->boundaries_special_treatment) {
int int_refinement = this->interpolate(point);
- Real refinement = this->interpolate(point)//int_refinement
+ Real refinement = int_refinement
* this->lattice->getRepresentativeConstant();
return distance_from_border >= refinement*.75; // heuristic value
} else {
return true;
}
}
}
template<Uint DIM>
void Domain<DIM>::compute_coupling_atoms(PointContainer<DIM> & interface,
PointContainer<DIM> & pad,
const Geometry<DIM> & atomistic) {
if (this->atoms.getSize() == 0) {
this->fill_atoms();
}
if (&atomistic == NULL) {
this->getMesh().compute_interface_atoms(interface, *this->refined);
} else {
this->getMesh().compute_interface_atoms(interface, atomistic);
}
// now, clean the pad atoms to get rid of the interface atoms
for (Uint interface_index = 0 ; interface_index < interface.getSize() ;
++interface_index) {
PointRef<DIM> interface_point = interface.getPoint(interface_index);
auto pad_index = this->pad_atoms_indices.begin();
auto pad_before_index = this->pad_atoms_indices.before_begin();
auto pad_end = this->pad_atoms_indices.end();
bool searching = true;
for (; pad_index != pad_end && searching; ++pad_index, ++pad_before_index) {
PointRef<DIM> pad_point = this->atoms.getPoint(*pad_index);
if (pad_point == interface_point) {
this->pad_atoms_indices.erase_after(pad_before_index);
searching=false;
}
}
}
auto pad_index = this->pad_atoms_indices.begin();
auto pad_end = this->pad_atoms_indices.end();
for(;pad_index != pad_end; ++pad_index) {
pad.addPoint(this->atoms.getPoint(*pad_index));
}
}
template class Domain<twoD>;
template class Domain<threeD>;
diff --git a/surface.org b/surface.org
index e1fc97e..8d21d89 100644
--- a/surface.org
+++ b/surface.org
@@ -1,10 +1,14 @@
* Surfaces Ids
when looping through surface elements
mesh. getuintdatapointer with tag_0
then push facets and tags together
(see mesh_io_msh)
* for reading (see single_edge_notched_beam.cc) decrement tag by one.
set tranction()
* in plugin:
- single_edge_notched_beam.cc:155
+ see single_edge_notched_beam.cc:155
+ add keywords in domain_akantu
+ - starting timestep
+ - surface_tag
+ - stress/force