diff --git a/src/common/point_container.cc b/src/common/point_container.cc index aa06499..4abb864 100644 --- a/src/common/point_container.cc +++ b/src/common/point_container.cc @@ -1,178 +1,185 @@ #include "point_container.hh" #include #include #include #include "geometry.hh" #include //template //using IntCoord = std::array; //template //using GridType = std::map, std::vector>; //template //using RevGridType = std::map>; template void PointContainer::filter(const Geometry & geom, Real tol) { PointContainer temp; for (Uint point = 0 ; point < this->getSize() ; ++point) { PointRef current_point = this->getPoint(point); if (geom.is_inside(current_point, tol)) { temp.addPoint(current_point); } } this->values = std::move(temp.values); this->modified = true; } // accumulates all the points of my box and all neighbouring boxes template inline void accumulate_test_points(const Uint & direction, const std::array & key, std::map, std::vector> & grid, std::array & off_set, std::vector & test_points) { for (int i = -1 ; i < 2 ; ++i) { off_set[direction] = key[direction] + i; if (direction < DIM -1) { accumulate_test_points(direction + 1, key, grid, off_set, test_points); } else { for (Uint point_id: grid[off_set]) { test_points.push_back(point_id); } } } } template void clean_grid(const Uint & index, std::map, std::vector> & grid, std::map> & reverse_grid) { std::array & box = reverse_grid[index]; auto pos = std::find(grid.at(box).begin(), grid.at(box).end(), index); grid.at(box).erase(pos); reverse_grid.erase(index); } template void PointContainer::clear() { this->threshold = -1; this->values.clear(); this->modified = true; } template Uint PointContainer::reportDuplicates(Real threshold) { this->setThreshold(threshold); return this->sniffDuplicates(); } template Uint PointContainer::cleanDuplicates(Real threshold) { this->setThreshold(threshold); return this->sniffDuplicates(); } template void PointContainer::createSearchGrid() { this->grid.clear(); this->reverse_grid.clear(); this->boxes.clear(); // sort all the points in a grid of size threshold std::array tmp_point; for (Uint point_id = 0 ; point_id < this->getSize() ; ++point_id) { for (Uint dim = 0 ; dim < DIM ; ++dim) { tmp_point[dim] = std::floor(this->getConstPoint(point_id)[dim]/this->threshold); } this->grid[tmp_point].push_back(point_id); this->reverse_grid[point_id] = tmp_point; } for (auto item: grid) { this->boxes.emplace_back(item.first); } - } template bool PointContainer::containsExacly(const PointRef & testpoint){ if (this->threshold <= 0) { this->setThreshold(1); } std::array tmp_point; for (Uint dim = 0 ; dim < DIM ; ++dim) { tmp_point[dim] = std::floor(testpoint[dim]/this->threshold); } for (auto point_id: this->grid[tmp_point]) { if (this->getConstPoint(point_id) == testpoint) { return true; } } return false; } template void PointContainer::setThreshold(Real thresh) { if (thresh <= 0.) { throw std::runtime_error("The threshold needs to be a positive finite real scalar"); } if ((thresh != this->threshold) || this->modified) { this->threshold = thresh; this->createSearchGrid(); } } template template Uint PointContainer::sniffDuplicates(){ Real tol2 = this->threshold*this->threshold; std::array tmp_point; // get potential neighbours for each point in the container + std::set nodes_to_delete; Uint counter = 0; for (auto key: this->boxes) { std::vector test_points; accumulate_test_points(0, key, this->grid, tmp_point, test_points); - for (Uint me : this->grid[key]) { - for (Uint other:test_points) { - if (me != other) { - Real r_square = distance2(this->getConstPoint(me), - this->getConstPoint(other)); - if (r_square < tol2) { - //There's a violation! kill it - counter++; - clean_grid(other, this->grid, this->reverse_grid); - if (do_delete) { - std::cout << "Deleting point " << other << " at position" << this->getConstPoint(other) << std::endl; - this->deletePoint(other); - this->threshold = -1.; //grid is invalidated - } else { - std::cout << "Duplicates: " << other << " at "; - this->getConstPoint(other).print_coords(std::cout); - std::cout << std::endl; - std::cout << " and: " << me << " at "; - this->getConstPoint(me).print_coords(std::cout); - std::cout << std::endl; - std::cout << " distance: " << std::sqrt(r_square) << std::endl; + std::vector & local_ids = this->grid[key]; + for (Uint me : local_ids) { + if (nodes_to_delete.find(me) == nodes_to_delete.end()) { + for (Uint other:test_points) { + if ((me != other) && + (nodes_to_delete.find(other) == nodes_to_delete.end())) { + Real r_square = distance2(this->getConstPoint(me), + this->getConstPoint(other)); + if (r_square < tol2) { + //There's a violation! kill it + counter++; + clean_grid(other, this->grid, this->reverse_grid); + if (do_delete) { + std::cout << "Deleting point " << other << " at position" << this->getConstPoint(other) << std::endl; + nodes_to_delete.insert(other); + } else { + std::cout << "Duplicates: " << other << " at "; + this->getConstPoint(other).print_coords(std::cout); + std::cout << std::endl; + std::cout << " and: " << me << " at "; + this->getConstPoint(me).print_coords(std::cout); + std::cout << std::endl; + std::cout << " distance: " << std::sqrt(r_square) << std::endl; + } } } } } } } if (do_delete) { if (counter > 0 ) { this->modified = true; + for (auto it = nodes_to_delete.rbegin(); it != nodes_to_delete.rend(); + ++it) { + this->deletePoint(*it); + } } } return counter; } template class PointContainer; template class PointContainer; diff --git a/src/common/point_container.hh b/src/common/point_container.hh index 39a54e0..66dc134 100644 --- a/src/common/point_container.hh +++ b/src/common/point_container.hh @@ -1,472 +1,479 @@ /** * @file point_container.hh * @author Till Junge * @date Fri Apr 13 11:10:22 2012 * * @brief * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #ifndef __CADD_MESH_POINT_CONTAINER_HH__ #define __CADD_MESH_POINT_CONTAINER_HH__ #include #include "string" #include "sstream" #include "fstream" #include "common.hh" #include "dumper_paraview.hh" #include "dumper_lammps.hh" #include #include #include +#include + template class Geometry; /*! \class PointRef simple representation of a point.*/ template class PointRef { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /*! Constructor * \param _vals array of which the point is part * \param _index position of the point in the array. E.g., for 3d points, the n-th point would have index n-1, and \a not 3*(n-1)*/ PointRef(const PointRef& other) :is_light(other.is_light) { if (this->is_light) { this->vals = other.vals; this->index = other.index; } else { this->vals = new Real [DIM]; std::copy(other.vals, other.vals + DIM, this->vals); this->index = 0; } }; PointRef & operator=(const PointRef& other) { if (this != &other) { if (this->is_light) { this->vals = new Real [DIM]; } std::copy(other.vals, other.vals + DIM, this->vals); this->is_light = false; this->index = 0; } return *this; } PointRef operator-(const PointRef& other) { PointRef tmp; std::copy(this->vals, this->vals + DIM, tmp.vals); for (Uint i = 0 ; i < DIM ; ++i) { tmp.vals[i] -= other[i]; } return tmp; } PointRef operator+(const PointRef& other) { PointRef tmp; std::copy(this->vals, this->vals + DIM, tmp.vals); for (Uint i = 0 ; i < DIM ; ++i) { tmp.vals[i] += other[i]; } return tmp; } PointRef(Real * _vals=NULL, Uint _index=0) : index(_index) { if (_vals == NULL) { this->is_light = false; this->vals = new Real[DIM]; } else { this->is_light = true; this->vals = _vals + _index * DIM; } } PointRef(std::initializer_list inlist) :vals(new Real[DIM]), index(0), is_light(false) { if (inlist.size() != DIM) { throw std::invalid_argument("Size of initializer list must be DIM"); } std::copy(inlist.begin(), inlist.end(), this->vals); } PointRef(std::vector & coords) :index(0), is_light(false) { this->vals = new Real[DIM]; for (Uint i = 0; i < DIM ; ++i) { this->vals[i] = coords.at(i); } } bool operator==(const PointRef & other) const { bool equal = true; for (Uint i = 0 ; i < DIM ; ++i) { equal &= (this->getComponent(i) == other.getComponent(i)); } return equal; } virtual ~PointRef(){ if (!this->is_light) delete[] this->vals;}; Real norm() {return vector_norm(this->vals);} /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void make_me_independent() { if (this-> is_light) { Real new_vals[DIM]; for (Uint i = 0 ; i < DIM ; ++i) { new_vals[i] = this->vals[i]; } this->vals = new_vals; this->is_light = false; } } void print_coords(std::ostream & stream) const { std::ios_base::fmtflags current_flags = stream.flags(); stream << std::scientific << "("; for (Uint i = 0; i < DIM -1; ++i) { stream << this->vals[i] << ", "; } stream << this->vals[DIM -1] << ")"; stream.setf(current_flags); } /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const { indent_space(stream, indent); stream << "Point "; stream.width(6); stream << this->index << ": "; this->print_coords(stream); } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /*! self-explanatory */ inline Real & operator[](Uint dir) {return this->vals[dir];} inline const Real & operator[](Uint dir) const {return this->vals[dir];} inline const Real & getComponent(Uint dir) const {return this->vals[dir];} inline Real * getArray() {return this->vals;} inline const Real * getArray() const {return this->vals;} inline std::vector getVector() const { return std::vector(this->vals, this->vals+DIM); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: Real * vals; Uint index; bool is_light; }; template inline Real distance2(const PointRef & pointA, const PointRef & pointB) { Real dist2 = 0; for (Uint i = 0; i < DIM ; ++i) { Real manhattan = pointA.getComponent(i) - pointB.getComponent(i); dist2 += manhattan * manhattan ; } return dist2; } template inline Real distance(const PointRef & pointA, const PointRef & pointB) { return std::sqrt(distance2(pointA, pointB)); } /*! \class PointContainer container for vectors */ template class PointContainer { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /// Constructor \param name by default "" PointContainer(std::string _name = ""):name(_name), threshold(-1.), modified(true){}; PointContainer(const PointContainer & other) :name(other.name), values(other.values), threshold(-1.), modified(true) {}; PointContainer(const std::vector & coords, const std::string & _name = "") :name(_name), values(coords), threshold(-1.), modified(true){}; virtual ~PointContainer(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: PointContainer & operator=(const PointContainer & other) { this->name = other.name; this->values = other.values; this->modified = true; return *this; } /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const { indent_space(stream, indent); stream << "PointContainer '" << this->name << "' contains " << this->getSize() << " points:" << std::endl; Uint size = this->getSize(); Uint max_size = 20; Uint i; for (i = 0; i < size && i < max_size; ++i) { PointRef printref = this->getConstPoint(i); printref.printself(stream, indent + indent_incr); stream << std::endl; } if (size > max_size) { indent_space(stream, indent + indent_incr); stream << "Point ..." << std::endl; } } void computeBounds(Real * bounds) const { for (Uint dir = 0 ; dir < DIM ; ++dir) { bounds[2 * dir + 0] = std::numeric_limits::max(); bounds[2 * dir + 1] = -std::numeric_limits::max(); } for (Uint point = 0 ; point < this->values.size() ; point+=DIM) { for (Uint dir = 0 ; dir < DIM ; ++dir) { keep_min(bounds[2*dir + 0], this->values.at(point+dir)); keep_max(bounds[2*dir + 1], this->values.at(point+dir)); } } } std::vector computeBounds() const { std::vector bounds; bounds.resize(2*DIM); this->computeBounds(&bounds[0]); return bounds; } void filter(const Geometry & geom, Real tol); void dump_paraview(std::string filename="diagnostic_pointcontainer") const { iohelper::DumperParaview dumper; dumper.setMode(iohelper::TEXT); dumper.setPoints(const_cast(&this->values[0]), DIM, this->values.size()/DIM, filename); dumper.setConnectivity(NULL, iohelper::POINT_SET, 0, iohelper::C_MODE); dumper.setPrefix("./"); dumper.init(); try { dumper.dump(); } catch (iohelper::IOHelperException & err) { std::cout << err.what() << std::endl; throw(err); } } void dump_lammps(std::string filename, std::vector bounds) const{ this->dump_lammps(filename, &bounds[0]); } void dump_lammps(std::string filename="diagnostic_pointcontainer", Real * set_bounds = NULL) const { Real bounds[2*DIM]; if (set_bounds == NULL) { this->computeBounds(bounds); } else { for (Uint i = 0 ; i < DIM*2 ; ++i) { bounds[i] = set_bounds[i]; } } iohelper::DumperLammps dumper(bounds); dumper.setPoints(const_cast(&this->values[0]), DIM, this->values.size()/DIM, filename); dumper.setPrefix("./"); dumper.setConnectivity(NULL, iohelper::POINT_SET, 0, iohelper::C_MODE); dumper.init(); dumper.dump(); } void dump_text(std::string filename="diagnostic_pointcontainer", std::string commentary="") const { std::ofstream ofile(filename.c_str()); if (ofile.good()) { ofile << "# file: " << filename <0) { ofile << "# " << commentary << std::endl; } } ofile << std::setprecision(16) << std::scientific << std::showpos; Uint current_point = 0; for (Uint i = 0 ; i < this->getSize(); ++i) { for (Uint j = 0 ; j < DIM - 1 ; ++j, ++current_point) { ofile << this->values.at(current_point) << "\t"; } ofile << this->values.at(current_point++); ofile << std::endl; } ofile.close(); } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /*! \param coords array of coordinates to add */ void addPoint(const Real * coords) { for (Uint i = 0; i < DIM; ++i) { this->values.push_back(coords[i]); } this->modified = true; } void addPoint(std::vector & coords) throw (std::string) { if (coords.size() < DIM) { std::stringstream error_msg; error_msg << "Size of coords should be " << DIM << " but it's " << coords.size(); throw error_msg.str(); } this->addPoint(&coords[0]); } inline void addPoint(const PointRef & point) { this->addPoint(&point.getComponent(0));} /// returns the number of points in the container inline Uint getSize() const {return this->values.size()/DIM;} /// returns the container name std::string getName() const {return this->name;} /// sets the container name void setName(std::string _name) {this->name = _name;} /*! returns a PointRef to the point at position \a index. Throws a string if * the index is out of bounds * \param index */ - PointRef getPoint(Uint index) throw (std::string) + PointRef getPoint(Uint index) throw (std::runtime_error) { - if (index > this->getSize()) { + if (index >= this->getSize()) { std::stringstream error_msg; error_msg << "index " << index << " is out of bounds. Size of vector " << this->name << " is only " << this->getSize(); - throw error_msg.str(); + throw std::runtime_error(error_msg.str()); } return PointRef(&this->values[0], index); } - const PointRef getConstPoint(Uint index) const throw (std::string) + const PointRef getConstPoint(Uint index) const throw (std::runtime_error) { - if (index > this->getSize()) { + if (index >= this->getSize()) { std::stringstream error_msg; error_msg << "index " << index << " is out of bounds. Size of vector " << this->name << " is only " << this->getSize(); - throw error_msg.str(); + throw std::runtime_error(error_msg.str()); } const PointRef ret_ref(const_cast(&values[0]), index); return ret_ref; } const PointRef operator[](Uint index) const throw (std::string) { return this->getConstPoint(index); } + PointRef operator[](Uint index) throw (std::string) { + return this->getConstPoint(index); + } /*! returnd the vector of raw data. You'll probably never use this.*/ const std::vector & getVector() const {return this-> values;} + const Real * getArray() const {return &this->values[0];} + /*! only use this if you know what you're doing. This can fuck up values and thus render the container useless */ std::vector & getNonConstVector() {return this->values;} void deletePoint(Uint index) { Uint nb_points = this->getSize(); for (Uint i = 0 ; i < DIM ; ++i) { this->values[DIM*index+i] = this->values[DIM*(nb_points-1) + i]; } for (Uint i = 0 ; i < DIM ; ++i) { this->values.pop_back(); } this->modified = true; } void clear(); Uint reportDuplicates(Real threshold); Uint cleanDuplicates(Real threshold); // check whether exactly this point (no rounding error) is in the container bool containsExacly(const PointRef &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: template Uint sniffDuplicates(); void createSearchGrid(); void setThreshold(Real thresh); std::string name; std::vector values; // search grid related members Real threshold; std::map, std::vector > grid; std::map > reverse_grid; std::vector > boxes; bool modified; //grid needs to be recomputed if modifications happened }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator <<(std::ostream & stream, const PointContainer & _this) { _this.printself(stream); return stream; } /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator <<(std::ostream & stream, const PointRef & _this) { _this.printself(stream); return stream; } #endif /* __CADD_MESH_POINT_CONTAINER_HH__ */ diff --git a/src/domain/CMakeLists.txt b/src/domain/CMakeLists.txt index 89c492a..2821eff 100644 --- a/src/domain/CMakeLists.txt +++ b/src/domain/CMakeLists.txt @@ -1,18 +1,19 @@ set(DOMAIN_SRC ${CMAKE_CURRENT_SOURCE_DIR}/tree.cc ${CMAKE_CURRENT_SOURCE_DIR}/cadd_mesh.cc ${CMAKE_CURRENT_SOURCE_DIR}/domain.cc ${CMAKE_CURRENT_SOURCE_DIR}/tet_mesh_quality.cc - ${CMAKE_CURRENT_SOURCE_DIR}/triangulation_quality.cc) + ${CMAKE_CURRENT_SOURCE_DIR}/triangulation_quality.cc + ${CMAKE_CURRENT_SOURCE_DIR}/mini_mesh.cc) set(DOMAIN_IGNORE_WARNINGS ${CMAKE_CURRENT_SOURCE_DIR}/tet_mesh_quality.cc ${CMAKE_CURRENT_SOURCE_DIR}/triangulation_quality.cc PARENT_SCOPE) set(CADD_MESH_SRC ${CADD_MESH_SRC} ${DOMAIN_SRC} PARENT_SCOPE) #add_library(domain ${DOMAIN_SRC}) #target_link_libraries(domain geometry common lattice akantu iohelper) diff --git a/src/domain/cadd_mesh.cc b/src/domain/cadd_mesh.cc index 29daf80..e9fd546 100644 --- a/src/domain/cadd_mesh.cc +++ b/src/domain/cadd_mesh.cc @@ -1,916 +1,1085 @@ extern "C" { #include #include #include #include #include /* for qh_vertexneighbors() */ #include /* for qh_facetarea(facet)*/ } #include "cadd_mesh.hh" #include #include #include #include #include #include #include #include "dumper_paraview.hh" // akantu includes for mesh io #include #include #include #include -// external routines for quality measures -#include "tet_mesh_quality.hh" -#include "triangulation_quality.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 Mesh::Mesh(const std::vector & points, Real tol, Geometry * 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 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 Mesh::~Mesh() { if (this->permutation != NULL) { delete this->permutation; delete this->remutation; } if (this->aka_mesh != NULL) { delete this->aka_mesh; } } template void Mesh::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 -inline Real computeRadiusRatio(const Real * points, const int * element, Uint nb_points) { - Uint nb_elements = 1; - Real q_min= std::numeric_limits::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 inline bool degenerate_finder(const std::vector & points, const std::vector & 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; const Real * points_ptr = &points.front(); const int * element_ptr = &element.front(); Real ratio = computeRadiusRatio(points_ptr, element_ptr, nb_nodes); // smaller ratio than criterion means element is good return (ratio < criterion) && (ratio > 0); } template inline void center_of_mass(const std::vector & points, const std::vector & element, PointRef & 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 inline void shear_points(std::vector & new_coords, const std::vector & 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 inline void smart_joggle(std::vector & new_coords, const std::vector & 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(new_coords, ref_coords, shear_dim, ref_dim, (ref_dim + shear_dim) * 1e-1); } } } template void Mesh::delaunay(const Real & degeneration_criterion, const std::string & flags){ // make a copy of the nodes before they get jiggled std::vector copy_nodes = this->nodes.getVector(); smart_joggle(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, ©_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 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 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(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 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 void Mesh::cleanup(const Geometry & atomistic_domain) { Uint nb_points = this->nodes.getSize(); Uint nb_elems = this->connectivity.size()/(DIM+1); std::vector new_points; std::vector new_connectivity; this->permutation = new std::vector(nb_points, -1); this->remutation = new std::vector(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(new_points); this->connectivity = new_connectivity; this->cleaned = true; } template void Mesh::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 Real Mesh::interpolate(const Real * coords, const std::vector & nodal_values) const { Real tol = 1e-6; 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 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::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::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(&(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 void Mesh::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 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 Real signOfJacobian(const Uint * elem, const std::vector & 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 void Mesh::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 & 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 * 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(this->connectivity.at((DIM+1)*elem_id + i)); } if (signOfJacobian(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 void Mesh::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 void Mesh::tagBoundaryPlane(Uint dir, Real val) { if (!this->aka_mesh_exists) { this->create_aka_mesh(); } akantu::MeshUtils::buildFacets(*this->aka_mesh); akantu::Array & nodes = this->aka_mesh->getNodes(); // find the elements of one dimension lower than the problem, they're the // facets typedef akantu::Array 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(*surface_type, "tag_0"); 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 inline bool compare_points(const PointRef & point1, const PointRef & point2) { bool retval = true; for (Uint i = 0; i < DIM ; ++i) { retval &= (fabs(1 - point1.getComponent(i)/point2.getComponent(i)) < std::numeric_limits::epsilon()); } return retval; } template void Mesh::compute_interface_atoms(PointContainer & container, const Geometry & atomistic, Real tol) { if (!this->aka_mesh_exists) { this->create_aka_mesh(); } akantu::MeshUtils::buildFacets(*this->aka_mesh); akantu::Array & nodes = this->aka_mesh->getNodes(); + // std::set interface_point_ids; + // // find the elements of one dimension lower than the problem, they're the + // // facets + // typedef akantu::Array 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) { + // PointRef point(&nodes((*element)(node))); + // if (atomistic.is_inside(point, tol)) { + // interface_point_ids.insert((*element)(node)); + // } + // } + // } + // } + // //the set is now filled with the indices of all interface atoms + // typedef std::set::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 point(&nodes(*point_id)); + // container.addPoint(point); + // } + + // get the set of all surface node ids std::set interface_point_ids; + this->surface_nodes(interface_point_ids); + typedef std::set::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 point(&nodes(*point_id)); + if (atomistic.is_inside(point, tol)) { + container.addPoint(point); + } + } +} + +/* -------------------------------------------------------------------------- */ +template +void Mesh::surface_nodes(std::set & 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 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) { - PointRef point(&nodes((*element)(node))); - if (atomistic.is_inside(point, tol)) { - interface_point_ids.insert((*element)(node)); - } + indices.insert((*element)(node)); } } } - //the set is now filled with the indices of all interface atoms - typedef std::set::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 point(&nodes(*point_id)); - container.addPoint(point); - } } template std::string Mesh::aka_mesh_name = "mesh 1"; template Uint Mesh::aka_mesh_counter = 1; double invert(double x){return 1./x;} template std::vector & Mesh::computeRadiusRatios(Uint index) { Uint nb_elements; int * tri_node; if (index == std::numeric_limits::max()) { - nb_elements = this->connectivity.size()/(DIM+1); + 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::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(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::iterator ratio = this->radius_ratios.begin(); std::vector::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 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 void Mesh::getDegenerateElements(const Real & quality_threshold, std::vector & indices) { if (!this->radius_ratio_computed) { this->computeRadiusRatios(); } greater_than_object comparator(quality_threshold); std::vector::const_iterator begin_ratio = this->radius_ratios.begin(); std::vector::iterator ratio = this->radius_ratios.begin(); std::vector::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 void Mesh::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 void Mesh::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(&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 void Mesh::splitDegenerates(const Real & radius_threshold, PointContainer & points) { typedef std::vector 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)); + 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 +void Mesh::get_neighbour_lists(std::vector> & 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::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." < +void Mesh::extract_neighbours(Uint elem_id, + const std::vector > & neighbours, + const std::set surf_ids, + const Geometry & atomic, + std::set & extracted_elem_ids, + std::set & free_nodes, + std::set & 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 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_degenerate = (degen_ids.find(neigh_id) != degen_ids.end()); + bool is_new = (extracted_elem_ids.find(neigh_id) == extracted_elem_ids.end()); + if (is_degenerate && is_new) { + this->extract_neighbours(neigh_id, neighbours, surf_ids, atomic, + extracted_elem_ids, free_nodes, degen_ids); + } else { + extracted_elem_ids.insert(neigh_id); + } + } + } +} + +/* -------------------------------------------------------------------------- */ +template +void Mesh::relaxDegenerates(const Real & radius_threshold, + const Geometry & atomic, + const Real & step_size, + const Real & min_step_size, + const Uint & max_iter, + bool blocked_nodes_lethal) { + std::set surf_nd_ids, degen_ids; + // get the set of surface nodes (they are always blocked) + this->surface_nodes(surf_nd_ids); + + typedef std::vector 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> neighbours; + this->get_neighbour_lists(neighbours); + + // loop through all the degenerates + while (!degen_ids.empty()) { + std::cout<< "There are " << degen_ids.size() + << " degenerates left to deal with" << std::endl; + std::set extracted_elem_ids; + std::set degenerate_neighbours; + std::set free_nodes; + + int current_elem_id = *degen_ids.begin(); + std::cout << "relaxing element " << current_elem_id + << " with quality threshold " << radius_threshold + << ". current quality is " + << computeRadiusRatio(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 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(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." < neighbourhood_connectivity; + for (const int & elem_id: extracted_elem_ids) { + neighbourhood_connectivity.push_back(this->getElem(elem_id)); + } + std::vector free_nodes_vec; + for (const int & node_id: free_nodes) { + free_nodes_vec.push_back(node_id); + } + MiniMesh neighbourhood(neighbourhood_connectivity, + free_nodes_vec, + this->nodes); + bool verbose = true; + try { + neighbourhood.relax(step_size, max_iter, radius_threshold, min_step_size, + verbose); + } catch (typename MiniMesh::StopIterError & e) { + if (blocked_nodes_lethal) { + throw e; + } else { + std::cout << "WARNING! " << e.what() << std::endl; + } + } + } +} +/* -------------------------------------------------------------------------- */ template class Mesh<2>; template class Mesh<3>; diff --git a/src/domain/cadd_mesh.hh b/src/domain/cadd_mesh.hh index 377c9e0..e55c174 100644 --- a/src/domain/cadd_mesh.hh +++ b/src/domain/cadd_mesh.hh @@ -1,236 +1,264 @@ /** * @file mesh.hh * @author Till Junge * @date Fri Apr 20 14:52:26 2012 * * @brief * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #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 +#include /*! \class Mesh builds FEM mesh */ template 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 & points, Real tol, Geometry * 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, const std::string & flags = default_flags); /*! 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 & 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 & nodal_values) const; inline Real interpolate(const PointRef & point, const std::vector & 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) ; bool containsExactNode(const PointRef & point) { return this->nodes.containsExacly(point);} /*! 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 coordinate 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 \param tol is the tolerance for deciding whether an atom is within atomistic */ void compute_interface_atoms(PointContainer & container, const Geometry & atomistic, Real tol); + /*! returns a set containing the surface node ids + */ + void surface_nodes(std::set & indices); /*! conputes the radius ratio for every element */ std::vector & computeRadiusRatios(Uint index = std::numeric_limits::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 & 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 & points); + /*! go from degenerate to degenerate and try to relax them by moving the nodes + * around. This functon is supposed to be smart and does not move atoms or + * surface points*/ + void relaxDegenerates(const Real & radius_threshold, + const Geometry & atomic, + const Real & step_size, + const Real & min_step_size, + const Uint & max_iter, + bool blocked_nodes_lethal = true); + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /*! returns the vector of nodes*/ const PointContainer & getNodes() const {return this->nodes;} /*! returns the vector of connectivities*/ const std::vector & getConnectivity() const {return this->connectivity;} + const int * getElem(int id) const {return &this->connectivity.at((DIM+1)*id);} + Uint getNbElems() const {return this->connectivity.size()/(DIM+1);} 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(); + // for each node, compute a list of all elements it appears in + void get_neighbour_lists(std::vector > & neighbours); + + // find and add all neighbouring elements to degenerate element elem_id. + // if a neighbour is itself also degenerate, recursively calls itself until done + void extract_neighbours(Uint elem_id, + const std::vector > & neighbours, + const std::set surf_ids, + const Geometry & atomic, + std::set & extracted_elem_ids, + std::set & free_nodes, + std::set & degen_ids); // contains all the nodes PointContainer nodes; // element definitions std::vector connectivity; // transformation matrices stored by matrix to evaluate interpolated nodal // values std::vector reverse_matrices; std::vector * permutation, * remutation; std::vector renumerotation; std::vector 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; Real tol; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator <<(std::ostream & stream, const Mesh & _this) { _this.printself(stream); return stream; } #endif /* __CADD_MESHER_MESH_HH__ */ diff --git a/src/domain/mesh_helper.hh b/src/domain/mesh_helper.hh new file mode 100644 index 0000000..d31c2db --- /dev/null +++ b/src/domain/mesh_helper.hh @@ -0,0 +1,28 @@ +// external routines for quality measures +#include "tet_mesh_quality.hh" +#include "triangulation_quality.hh" + +template +inline Real computeRadiusRatio(const Real * points, const int * element, Uint nb_points) { + Uint nb_elements = 1; + Real q_min= std::numeric_limits::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; +} + diff --git a/src/domain/mini_mesh.cc b/src/domain/mini_mesh.cc new file mode 100644 index 0000000..e7cddff --- /dev/null +++ b/src/domain/mini_mesh.cc @@ -0,0 +1,118 @@ +#include "mini_mesh.hh" +#include "mesh_helper.hh" +#include "common.hh" +#include + + +/* -------------------------------------------------------------------------- */ +template +MiniMesh::MiniMesh(const std::vector & connectivity, + const std::vector & free_nodes, + PointContainer & nodes) + :connectivity(connectivity), free_nodes(free_nodes), nodes(nodes), + grad(DIM*free_nodes.size()) {} + +/* -------------------------------------------------------------------------- */ +template +Real MiniMesh::computeQuality() { + Real quality = 1; + for (auto elem: this->connectivity) { + keep_max(quality, + computeRadiusRatio(this->nodes.getArray(), elem, + this->nodes.getSize())); + } + return quality; +} + +/* -------------------------------------------------------------------------- */ +template +void MiniMesh::updateGradient(const Real & step_size) { + const Real ref_quality = this->computeQuality(); + Uint nb_free_nodes = this->free_nodes.size(); + for (Uint local_nd_id = 0 ; local_nd_id < nb_free_nodes ; ++local_nd_id) { + Uint node_id = this->free_nodes[local_nd_id]; + for (Uint i = 0 ; i < DIM ; ++i) { + this->nodes[node_id][i] += step_size; + const Real new_quality = this->computeQuality(); + this->nodes[node_id][i] -= step_size; + this->grad.at(DIM*local_nd_id+i) = (new_quality-ref_quality)/step_size; + } + } +} + +/* -------------------------------------------------------------------------- */ +template +void MiniMesh::normGradient(Real step_size) { + Real max_step = 0; + Uint nb_free_nodes = this->free_nodes.size(); + for (Uint step_id = 0 ; step_id < nb_free_nodes ; ++step_id) { + keep_max(max_step, vector_norm(&this->grad[DIM*step_id])); + } + Real factor = step_size/max_step; + for (Real & val: this->grad) { + val *= factor; + } +} + +/* -------------------------------------------------------------------------- */ +template +Uint MiniMesh::relax(Real initial_step_size, + Uint max_iter, Real threshold, + Real min_step_size, + bool verbose) { + Uint iter = 0; + Real quality = this->computeQuality(); + Real old_quality = quality + 1; + Uint nb_free_nodes = this->free_nodes.size(); + Uint active_max_iter = max_iter; + while (quality > threshold) { + if (verbose) { + std::cout << " At beginning of iteration " << iter+1 + << ", local quality is " << quality << std::endl; + } + if(++iter >= active_max_iter) { + throw StopIterError("reached max_iter without reaching threshold"); + } + this->updateGradient(initial_step_size); + this->normGradient(initial_step_size); + for (Uint local_nd_id = 0 ; local_nd_id < nb_free_nodes ; ++local_nd_id) { + Uint global_nd_id = this->free_nodes[local_nd_id]; + for (Uint i = 0 ; i < DIM ; ++i) { + this->nodes[global_nd_id][i] -= this->grad[DIM*local_nd_id+i]; + } + } + old_quality = quality; + quality = this->computeQuality(); + // make sure it was an improvement + if ((old_quality < quality) && (initial_step_size > min_step_size)) { + quality = old_quality; + // we have to roll back and reduce the step size + for (Uint local_nd_id = 0 ; local_nd_id < nb_free_nodes ; ++local_nd_id) { + Uint global_nd_id = this->free_nodes[local_nd_id]; + for (Uint i = 0 ; i < DIM ; ++i) { + this->nodes[global_nd_id][i] += this->grad[DIM*local_nd_id+i]; + } + } + initial_step_size *= .5; + keep_max(initial_step_size, min_step_size); + if (verbose) { + std::cout << "Reduced step size to " << initial_step_size << endl; + } + active_max_iter = max_iter + iter; + } + + } + if (verbose) { + std::cout << " Success: reached a quality of " << quality << " in " << iter + << " iteration"; + if (iter != 1){ + std::cout << "s"; + } + std::cout << std::endl; + } + return iter; +} + +/* -------------------------------------------------------------------------- */ +template class MiniMesh; +template class MiniMesh; diff --git a/src/domain/mini_mesh.hh b/src/domain/mini_mesh.hh new file mode 100644 index 0000000..2bd97ed --- /dev/null +++ b/src/domain/mini_mesh.hh @@ -0,0 +1,49 @@ + +#ifndef _MINI_MESH_H_ +#define _MINI_MESH_H_ + +#include "common.hh" +#include "point_container.hh" +#include + + +/*! This class allows to create miniature meshes based on a connectivity and a + * point container for quality calculations and similar + */ +template +class MiniMesh { +public: + MiniMesh(const std::vector & connectivity, + const std::vector & free_nodes, + PointContainer & nodes); + virtual ~MiniMesh(){}; + class StopIterError: public std::runtime_error { + public: + StopIterError(std::string str):std::runtime_error(str){}; + virtual ~StopIterError(){}; + }; + + Uint relax(Real imposed_step_size, + Uint max_iter, + Real threshold, + Real min_step_size, + bool verbose = false); + +private: + /// returns the maximum quality (worst element) in the minimesh + Real computeQuality(); + + /// computes numerically the gradient of the quality + void updateGradient(const Real & step_size); + + /// normalize the gradient so that the largest point move is step_size + void normGradient(Real step_size); + + const std::vector & connectivity; + const std::vector & free_nodes; + PointContainer & nodes; + + std::vector grad; +}; + +#endif /* _MINI_MESH_H_ */ diff --git a/src/geometry/cylinder.cc b/src/geometry/cylinder.cc index e963010..262fca3 100644 --- a/src/geometry/cylinder.cc +++ b/src/geometry/cylinder.cc @@ -1,385 +1,388 @@ #include "cylinder.hh" #include #include #include #include #include /* -------------------------------------------------------------------------- */ class CylinderException: public std::exception { public: CylinderException(std::stringstream & _what): what_stream(_what.str()) {}; CylinderException(std::string & _what): what_stream(_what) {}; ~CylinderException() throw() {}; virtual const char* what() const throw() { return this->what_stream.c_str(); } private: std::string what_stream; }; /* -------------------------------------------------------------------------- */ template Cylinder::Cylinder(const std::vector & centre, const Real & radius_, const std::vector & axis_, std::string name) :Geometry(name) { std::stringstream error; if (DIM < 3) { error << "This constructor only makes sense for more than 2 dimensions."; throw CylinderException(error); } Uint size_centre = centre.size(); if (size_centre != DIM) { error << "The length of centre should have been " << DIM << " but I got " << size_centre << "." ; throw CylinderException(error); } Uint size_axis = axis_.size(); if (size_axis != DIM) { error << "The length of axis should have been " << DIM << " but I got " << size_centre << "." ; throw CylinderException(error); } this->length = vector_norm(&axis_[0]); for (Uint dir = 0 ; dir < DIM ; ++dir) { this->base[dir] = centre[dir]; this->axis[dir] = axis_[dir]/this->length; } this->radius = radius_; } /* -------------------------------------------------------------------------- */ template Cylinder::Cylinder(const std::vector & centre, const Real & radius_, std::string name) :Geometry(name) { std::stringstream error; if (DIM > 2) { error << "This constructor only makes sense for less than 3 dimesions."; throw CylinderException(error); } Uint size_centre = centre.size(); if (size_centre != DIM) { error << "The length of centre should have been " << DIM << " but I got " << size_centre << "." ; throw CylinderException(error); } for (Uint dir = 0 ; dir < DIM ; ++dir) { this->base[dir] = centre[dir]; this->axis[dir] = 0; } this->radius = radius_; this->length = 1.;//vector_norm(this->axis); } /* -------------------------------------------------------------------------- */ template Cylinder::Cylinder(const PointRef & centre, const Real & radius, const Real * axis_, std::string name) :Geometry(name) { this->base = centre; this->length = vector_norm(axis_); for (Uint dir = 0 ; dir < DIM ; ++dir) { this->axis[dir] = axis_[dir]/this->length; } this->radius = radius; } /* -------------------------------------------------------------------------- */ template Cylinder::Cylinder(const Cylinder & other) { vector_set(this->base.getArray(), other.base.getArray()); vector_set(this->axis, other.axis); this->length = other.length; this->radius = other.radius; } /* -------------------------------------------------------------------------- */ template Cylinder::~Cylinder(){}; /* -------------------------------------------------------------------------- */ template void Cylinder::printself(std::ostream & stream, int indent) const { indent_space(stream, indent); stream << "Cylinder " << this->get_name() << ": " << std::endl; indent_space(stream, indent); stream << "base : " << this->base << std::endl; indent_space(stream, indent); stream << "axis : " << PointRef(const_cast(this->axis)) << std::endl; indent_space(stream, indent); stream << "radius : " << this->radius; } /* -------------------------------------------------------------------------- */ template void Cylinder::shift(const PointRef & offset) { for (Uint i = 0 ; i < DIM ; ++i) { this->base[i] += offset[i]; } } /* -------------------------------------------------------------------------- */ template InsideObject Cylinder::from_border(const Real * point, Real tol) const { PointRef relpos = PointRef(const_cast(point)) - this->base; Real violation = norm_cross_prod(relpos.getArray(), axis) - this->radius; if (DIM == threeD) { Real dot_p = dot_prod(relpos.getArray(), axis); keep_max(violation, -dot_p); keep_max(violation, dot_p - this->length); } return {violation, violation <= 0}; } /* -------------------------------------------------------------------------- */ template Real Cylinder::min_in_direction(const Real * dir, const Real & unit) const { Real norm = -vector_norm(dir); Real dist = this->in_direction(dir, norm); return dist/unit; } /* -------------------------------------------------------------------------- */ template Real Cylinder::max_in_direction(const Real * dir, const Real & unit) const { Real norm = vector_norm(dir); Real dist = this->in_direction(dir, norm); return dist/unit; } /* -------------------------------------------------------------------------- */ template Real Cylinder::size_in_direction(const Real * dir, const Real & unit) const { Real norm = vector_norm(dir); Real dist = this->in_direction(dir, norm) - this->in_direction(dir, -norm); return dist/unit; } /* -------------------------------------------------------------------------- */ template inline void circle_filler(const Real r1_unit[DIM], const Real r2_unit[DIM], const Real centre[DIM], const Real & radius, const Mesh & auxiliary, const std::vector & refinement, const Real step_size[DIM], const Real & repr_const, const bool periodicity[DIM], PointContainer & points) { Real alpha = 0.; while (alpha < 2*pi) { Real periphery[DIM] = {0}; vector_incr(periphery, centre); vector_incr(periphery, r1_unit, cos(alpha) * radius); vector_incr(periphery, r2_unit, sin(alpha) * radius); PointRef tmp(periphery); add_if_ok(tmp, auxiliary, refinement, repr_const, points); // a bit of trig Real delta_alpha = 2 * acos(sqrt(1-square(repr_const)/ square(radius))); alpha += delta_alpha; } } /* -------------------------------------------------------------------------- */ template void Cylinder::generate_surface_points(const Mesh & auxiliary, const std::vector & refinement, const Real step_size[DIM], const Real & repr_const, const bool periodicity[DIM], PointContainer & points) const { this->periodicity_check(periodicity); - Real r1_unit[DIM] = {1., 0.}; - Real r2_unit[DIM] = {0., 1.}; + Real r1_unit[DIM] = { 0.}; + r1_unit[0] = 1.; + Real r2_unit[DIM] = { 0.}; + r2_unit[1] = 1.; if (DIM == twoD) { circle_filler(r1_unit, r2_unit, this->base.getArray(), this->radius, auxiliary, refinement, step_size, repr_const, periodicity, points); } else if (DIM == threeD) { // compute a vector of length R in the plane of the base - Real helper[DIM]; + Real helper[DIM] = {0.}; for (Uint i = 0 ; i < DIM ; ++i) { helper[(i+1)%DIM] = this->axis[i]; } cross_prod(this->axis, helper, r1_unit); cross_prod(this->axis, r1_unit, r2_unit); - Real top[DIM], bottom[DIM] = {0}; + Real top[DIM] = {0.}, bottom[DIM] = {0.}; vector_incr(top, this->axis, this->length); vector_incr(top, this->base.getArray()); //start by filling the top and bottom of the cylinder for (Real fill_radius = this->radius; fill_radius > 0; fill_radius -= repr_const) { circle_filler(r1_unit, r2_unit, this->base.getArray(), fill_radius, auxiliary, refinement, step_size, repr_const, periodicity, points); circle_filler(r1_unit, r2_unit, top, fill_radius, auxiliary, refinement, step_size, repr_const, periodicity, points); fill_radius -= repr_const; } //fill the cylindrical surface vector_set(bottom, this->base.getArray()); for (Real offset = 0 ; offset < this->length/2. ; offset += repr_const) { //set the offsets for the filler vector_incr(bottom, this->axis, repr_const); vector_incr(top, this->axis, -repr_const); // fill the ring circle_filler(r1_unit, r2_unit, bottom, this->radius, auxiliary, refinement, step_size, repr_const, periodicity, points); circle_filler(r1_unit, r2_unit, top, this->radius, auxiliary, refinement, step_size, repr_const, periodicity, points); } } else { std::stringstream error; error << "Not defined for DIM = " << DIM << "."; throw CylinderException(error); } + points.cleanDuplicates(1); } /* -------------------------------------------------------------------------- */ template void Cylinder::complement_periodically(const Mesh & auxiliary, const std::vector & refinement, PointContainer & points, int direction, Uint dim) const { bool periodicity[DIM] = {0}; periodicity[dim] = true; this->periodicity_check(periodicity); Real offset[DIM] = {0}; /* positive direction means that bottom points are replicated at the top and * vice-versa */ offset[dim] = std::copysign(this->length, direction); Real reference_point[DIM]; vector_set(reference_point, this->base.getArray()); if (direction < 0) { vector_incr(reference_point, offset, 1.); } vector_incr(offset, reference_point); PointRef point; Real tol = std::max(fabs(reference_point[dim]), 1.); tol = 2000*tol*std::numeric_limits::epsilon(); Uint nb_points = points.getSize(); for (Uint i = 0; i < nb_points ; ++i) { point = points.getPoint(i); if (fabs(point.getComponent(dim) - reference_point[dim]) < tol) { point[dim] = offset[dim]; add_periodic_complement(auxiliary, refinement, point, points); } } } /* -------------------------------------------------------------------------- */ template inline Real circle_extreme(const Real * centre, const Real * axis_unit, const Real * dir_norm, const Real & radius) { std::stringstream error; error << "circle_extreme is not implemented for DIM = " << DIM << "."; throw CylinderException(error); } /* -------------------------------------------------------------------------- */ template <> inline Real circle_extreme(const Real * centre, const Real * axis_unit, const Real * dir_norm, const Real & radius) { const Uint DIM = threeD; Real helper[DIM], r1_unit[DIM], r2_unit[DIM]; for (Uint i = 0 ; i < DIM ; ++i) { helper[(i+1)%DIM] = axis_unit[i]; } cross_prod(axis_unit, helper, r1_unit); cross_prod(axis_unit, r1_unit, r2_unit); Real cos1 = dot_prod(dir_norm, r1_unit); Real cos2 = dot_prod(dir_norm, r2_unit); Real extreme[DIM]; for (Uint i = 0 ; i < DIM ; ++i) { extreme[i] = centre[i] + radius * (r1_unit[i] * cos1 + r2_unit[i] * cos2); } return dot_prod(extreme, dir_norm); } /* -------------------------------------------------------------------------- */ template <> inline Real circle_extreme(const Real * centre, const Real * axis_unit, const Real * dir_norm, const Real & radius) { const Uint DIM = twoD; Real helper[DIM]; vector_set(helper, centre); vector_incr(helper, dir_norm, radius); return dot_prod(helper, dir_norm); } /* -------------------------------------------------------------------------- */ template Real Cylinder::in_direction(const Real * dir, const Real & norm) const { Real dir_norm[DIM]; vector_set(dir_norm, dir); vector_scale(dir_norm, 1./norm); // check base Real dist_base = circle_extreme(this->base.getArray(), this->axis, dir_norm, this->radius); if (DIM == twoD) { return dist_base; } // check top Real top[DIM]; vector_set(top, this->base.getArray()); vector_incr(top, this->axis, this->length); Real dist_top = circle_extreme(top, this->axis, dir_norm, this->radius); if (norm*dist_base > norm*dist_top) { return dist_base; } return dist_top; } /* -------------------------------------------------------------------------- */ template void Cylinder::periodicity_check(const bool periodicity[DIM]) const { Real err2 = 0; bool periodicity_exists = false; for (Uint i = 0 ; i < DIM ; ++i) { err2 += this->axis[i] - static_cast(periodicity[i]); periodicity_exists |= periodicity[i]; } if (periodicity_exists && (err2 > 8*std::numeric_limits::epsilon())) { std::stringstream error; error << "Periodicity (if any) MUST be aligned with the cylinder axis. " << "This axis is "; print_array(this->axis, DIM, "axis", error); error << std::endl << " and the periodicity is "; print_array(periodicity, DIM, "periodicity", error); throw CylinderException(error); } } template class Cylinder; template class Cylinder; diff --git a/src/swig/CaddMesher.i b/src/swig/CaddMesher.i index e699316..b27ad03 100644 --- a/src/swig/CaddMesher.i +++ b/src/swig/CaddMesher.i @@ -1,169 +1,170 @@ %module CADDMesher #pragma SWIG nowarn=325,389 %{ #include "../src/common/cadd_mesher.hh" %} typedef double Real; typedef unsigned int Uint; %include "std_string.i" %include "std_vector.i" +%include "std_set.i" %include "typemaps.i" %define __STR__() std::string __str__() { std::ostringstream out; out << *$self; return out.str(); } %enddef namespace std { %template(vecReal) std::vector ; %template(vecInt) std::vector ; %template(vecBool) std::vector ; } %rename(__assign__) *::operator=; %apply Real &INPUT {Real & refinement}; %include "../src/domain/domain.hh" //-------------------------------------------- // Domain %extend Domain { __STR__()}; %ignore Domain::setPointRefinement(Real const *, Real &); %template(Domain2d) Domain; %extend Domain { __STR__()}; %ignore Domain::setPointRefinement(Real const *, Real &); %template(Domain3d) Domain; //-------------------------------------------- // Geometry %include "../src/geometry/geometry.hh" %extend Geometry { __STR__() }; %template(Geometry2d) Geometry; %extend Geometry { __STR__() }; %template(Geometry3d) Geometry; %include "../src/geometry/block.hh" %ignore Block::Block( Real const *, std::string); %template(Block2d) Block; %ignore Block::Block( Real const *, std::string); %template(Block3d) Block; %include "../src/geometry/cylinder.hh" %ignore Cylinder::Cylinder( Real const *, std::string); %template(Cylinder2d) Cylinder; %ignore Cylinder::Cylinder( Real const *, std::string); %template(Cylinder3d) Cylinder; //-------------------------------------------- // Lattice %include "../src/lattice/lattice.hh" %extend Lattice { __STR__() }; %template(Lattice2d) Lattice; %extend Lattice { __STR__() }; %template(Lattice3d) Lattice; %include "../src/lattice/fcc_lattice.hh" %ignore FccLattice::FccLattice(Real *, Real *); %ignore FccLattice::FccLattice(Real *); %extend FccLattice { __STR__() }; %template(FccLattice2d) FccLattice; %ignore FccLattice::FccLattice(Real *, Real *); %ignore FccLattice::FccLattice(Real *); %extend FccLattice { __STR__() }; %template(FccLattice3d) FccLattice; //-------------------------------------------- // PointContainer and PointRefs %include "../src/common/point_container.hh" %ignore PointContainer::addPoint(const Real*); %extend PointContainer { __STR__() }; %template(PointContainer2d) PointContainer; %ignore PointContainer::addPoint(const Real*); %extend PointContainer { __STR__() }; %template(PointContainer3d) PointContainer; %ignore PointRef::operator=; %ignore PointRef::PointRef(Real *, Uint); %ignore PointRef::PointRef(Real *); %extend PointRef { __STR__()}; %extend PointRef { inline void __setitem__(size_t i, const Real & val) throw(std::out_of_range) { if (i>=twoD || i<0) { throw std::out_of_range("out of bounds access"); } $self->getArray()[i] = val; } } %extend PointRef { inline double __getitem__(size_t i) throw(std::out_of_range) { if (i>=twoD || i<0) { throw std::out_of_range("out of bounds access"); } return $self->getArray()[i]; } } %template (PointRef2d) PointRef; %ignore PointRef::PointRef(Real *, Uint); %ignore PointRef::PointRef(Real *); %extend PointRef { __STR__()}; %extend PointRef { inline void __setitem__(size_t i, const Real & val) throw(std::out_of_range) { if (i>=threeD || i<0) { throw std::out_of_range("out of bounds access"); } $self->getArray()[i] = val; } } %extend PointRef { inline double __getitem__(size_t i) throw(std::out_of_range) { if (i>=threeD || i<0) { throw std::out_of_range("out of bounds access"); } return $self->getArray()[i]; } } %template (PointRef3d) PointRef; //-------------------------------------------- // Mesh %include "../src/domain/cadd_mesh.hh" %extend Mesh { __STR__()}; %template (Mesh2d) Mesh; %extend Mesh { __STR__()}; %template (Mesh3d) Mesh; //-------------------------------------------- // Tree %include "../src/domain/tree.hh" %extend Tree { __STR__()}; %ignore Tree::iterator; %ignore Tree::getIterator(); %template (Tree2d) Tree; %extend Tree { __STR__()}; %ignore Tree::iterator; %ignore Tree::getIterator(); %template (Tree3d) Tree; %include "../src/common/cadd_mesher.hh" diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5987170..83b51d8 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,54 +1,59 @@ #/** #* @file CMakeLists.txt #* @author Till Junge #* @date Thu Apr 12 15:05:28 2012 #* #* @brief #* #* @section LICENSE #* #* #* #*/ # #/* -------------------------------------------------------------------------- */ option(BLOCK_TEST "build the block test") if (BLOCK_TEST) register_test(block_test block_test.cc) endif(BLOCK_TEST) option(POINT_CONTAINER_TEST "build the point_container test") if (POINT_CONTAINER_TEST) register_test(point_container_test point_container_test.cc) endif(POINT_CONTAINER_TEST) option(FCC_LATTICE_TEST "build the fcc_lattice test") if (FCC_LATTICE_TEST) register_test(fcc_lattice_test fcc_lattice_test.cc) endif(FCC_LATTICE_TEST) option(MESH_TEST "build the mesh test") if (MESH_TEST) register_test(mesh_test mesh_test.cc) endif(MESH_TEST) option(MESH_TEST_2D "build the 2D mesh test") if (MESH_TEST_2D) register_test(mesh_test_2d mesh_test_2d.cc) endif(MESH_TEST_2D) option(DOMAIN_TEST "build the domain test") if (DOMAIN_TEST) register_test(domain_test domain_test.cc) endif(DOMAIN_TEST) option(DOMAIN_TEST_3D "build the 3D domain test") if (DOMAIN_TEST_3D) register_test(domain_test_3d domain_test_3d.cc) endif(DOMAIN_TEST_3D) option(ROTATION_TEST_3d "build the 3d rotation test") if (DOMAIN_TEST_3D) register_test(rotation_test_3d rotation_test_3d.cc) endif(DOMAIN_TEST_3D) + +option(MINI_MESH_TEST "Test relaxation of degenerates") +if (MINI_MESH_TEST) + register_test(mini_mesh_test mini_mesh_test.cc) +endif(MINI_MESH_TEST) diff --git a/test/mini_mesh_test.cc b/test/mini_mesh_test.cc new file mode 100644 index 0000000..b691073 --- /dev/null +++ b/test/mini_mesh_test.cc @@ -0,0 +1,61 @@ +#include +#include "common.hh" +#include "mini_mesh.hh" +#include "point_container.hh" +#include + + +const Uint dim = 2; +const Uint nb_elem = 8; +int conn[nb_elem][dim+1] = {{0, 5, 4}, + {0, 1, 5}, + {1, 6, 5}, + {1, 2, 6}, + {2, 4, 6}, + {2, 3, 4}, + {3, 0, 4}, + {4, 5, 6}}; + +void createNodes(PointContainer & nodes, std::vector & free_nodes) { + nodes.addPoint({0., 0. }); + nodes.addPoint({4., 0. }); + nodes.addPoint({4., 2. }); + nodes.addPoint({0., 2. }); + nodes.addPoint({1., 1.1}); + nodes.addPoint({2., 0.9}); + nodes.addPoint({3., 1.1}); + free_nodes.push_back(4); + free_nodes.push_back(5); + free_nodes.push_back(6); +} + +void createConnectivity(std::vector & connectivity) { + for (Uint i = 0 ; i < nb_elem ; ++i) { + connectivity.push_back(conn[i]); + } +} + +int main(int argc, char *argv[]) +{ + PointContainer nodes("Nodes"); + std::vector free_nodes; + createNodes(nodes, free_nodes); + + std::vector connectivity; + createConnectivity(connectivity); + + MiniMesh mesh(connectivity, free_nodes, nodes); + Uint iter = mesh.relax(.1, 10, 10, true); + std::cout<< nodes << std::endl; + + Uint expected_iter = 1; + if (iter == expected_iter){ + std::cout << "Converged in " << iter + << " steps. This is as it should be. SUCCESS." << std::endl; + return 0; + } else { + std::cout << "This should have converged in " << expected_iter << " steps " + << " but did so in " << iter << ". FAILURE." << std::endl; + return -1; + } +} diff --git a/test/point_container_test.cc b/test/point_container_test.cc index 58e9498..f793185 100644 --- a/test/point_container_test.cc +++ b/test/point_container_test.cc @@ -1,76 +1,93 @@ /** * @file point_container_test.cc * @author Till Junge * @date Fri Apr 13 13:00:14 2012 * * @brief * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #include #include "common.hh" #include "point_container.hh" int main(int argc, char *argv[]) { const Uint nb_dat = 24; Real data[nb_dat]; for (Uint i = 0; i < nb_dat; ++i) { data[i] = 12.2*i; } Uint dim = 2; PointRef test_ref(data, twoD); PointRef test_copy; test_copy = test_ref; std::cout << "test_ref = " << test_ref << ", test_copy = " << test_copy << std::endl; PointContainer<2> cont2("2D example"); for (Uint i = 0; i < nb_dat/dim; ++i) { cont2.addPoint(data + dim*i); } std::cout << "container size should be " << nb_dat/dim << " and it is " << cont2.getSize() << std::endl; if (nb_dat/dim != cont2.getSize()) { std::cout << "Test Failed" < cont3("3D example"); for (Uint i = 0; i < nb_dat/dim; ++i) { cont3.addPoint(data + dim*i); } std::cout << "container size should be " << nb_dat/dim << " and it is " << cont3.getSize() << std::endl; if (nb_dat/dim != cont3.getSize()) { std::cout << "Test Failed" < duplicate_container; + duplicate_container.addPoint({0, 0}); + duplicate_container.addPoint({2, 0}); + duplicate_container.addPoint({2, 2}); + duplicate_container.addPoint({0, 2}); + duplicate_container.addPoint({1, 1}); + duplicate_container.addPoint({1, 1}); + + duplicate_container.cleanDuplicates(1.); + + duplicate_container.clear(); + duplicate_container.addPoint({0, 0}); + duplicate_container.addPoint({2, 0}); + duplicate_container.addPoint({0, 2}); + duplicate_container.addPoint({1, 1}); + duplicate_container.addPoint({1, 1}); + duplicate_container.addPoint({2, 2}); + + duplicate_container.cleanDuplicates(1.); std::cout << "Test seems to have passed. Please check output." << std::endl; exit(0); }