diff --git a/src/domain/cadd_mesh.cc b/src/domain/cadd_mesh.cc index e9fd546..7c0bb0b 100644 --- a/src/domain/cadd_mesh.cc +++ b/src/domain/cadd_mesh.cc @@ -1,1085 +1,1090 @@ 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 #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 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]={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) { indices.insert((*element)(node)); } } } } 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->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}; 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) { + bool blocked_nodes_lethal, + bool verbose) { 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; + if (verbose) { + 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; + if (verbose) { + 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; + std::cout << "WARNING! Element " <; template class Mesh<3>; diff --git a/src/domain/cadd_mesh.hh b/src/domain/cadd_mesh.hh index e55c174..e47c8f8 100644 --- a/src/domain/cadd_mesh.hh +++ b/src/domain/cadd_mesh.hh @@ -1,264 +1,265 @@ /** * @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); + bool blocked_nodes_lethal = true, + bool verbose = false); /* ------------------------------------------------------------------------ */ /* 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/domain.cc b/src/domain/domain.cc index 8d0ca36..40ba574 100644 --- a/src/domain/domain.cc +++ b/src/domain/domain.cc @@ -1,263 +1,276 @@ #include "domain.hh" #include "tree.hh" #include #include #include #include template Domain::Domain(const Geometry & overall_, const Geometry & refined_, const Geometry & atomistic_, const Geometry & atomistic_skinned_, const Lattice & lattice_, bool boundaries_special_treatment_, const PointRef & centre, Real tolerance, const std::string & qh_flags_) :auxiliary_mesh(NULL), tree(NULL), main_mesh(NULL), filled(false), main_meshed(false), boundaries_special_treatment(boundaries_special_treatment_), mesh_quality(-std::numeric_limits::max()), tol(tolerance), qh_flags(qh_flags_) { 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::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(NULL, 2*pow2size, this->lattice, this); } else { this->tree = new Tree(centre.getArray(), 2*pow2size, this->lattice, this); } } template Domain::~Domain() { delete this->auxiliary_mesh; delete this->main_mesh; delete this->tree; } template void Domain::setPointRefinement(const Real * coords, Real & refinement ) { this->refinement_points.addPoint(coords); this->refinement.push_back(refinement); } template void Domain::fill_points(bool periodicity[DIM]) throw(std::string){ this->initialiseAuxiliaryMesh(); this->auxiliary_mesh->dump_paraview("test"); if (this->boundaries_special_treatment) { this->overall->generate_surface_points(*this->auxiliary_mesh, this->refinement, this->lattice->getSpatialConstants(), this->lattice->getRepresentativeConstant(), periodicity, this->points); } this->tree->fill(this->points); this->filled = true; this->points.reportDuplicates(1); } template void Domain::complement_periodically(int direction, Uint dim) throw(std::string){ if (!this->filled) { throw("Mesh has got to be filled before it can be complemented"); } this->overall->complement_periodically(*this->auxiliary_mesh, this->refinement, this->points, direction, dim); } template void Domain::printself(std::ostream & stream, int indent) const{ indent_space(stream, indent); stream << "PRINTSELF IS NOT YET IMOPLEMENTED FOR CLASS DOMAIN"; } template void Domain::initialiseAuxiliaryMesh() { if (this->auxiliary_mesh != NULL) { delete this->auxiliary_mesh; } this->auxiliary_mesh = new Mesh(this->refinement_points.getVector(), 0); this->auxiliary_mesh->delaunay(1000, this->qh_flags); this->auxiliary_mesh->initInterpolation(); } template Mesh & Domain::getMesh(Real quality) throw(std::string) { if ((!this->main_meshed) || (quality!=this->mesh_quality)) { this->main_mesh = new Mesh(this->points.getVector(), this->tol, this->atomistic_skinned); try { this->main_mesh->delaunay(quality, this->qh_flags); } 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; this->fill_atoms(); this->fill_coupling_atoms(); this->mesh_quality = quality; } return *this->main_mesh; } template void Domain::splitDegenerates(const Real & radius_threshold) { this->main_mesh->splitDegenerates(radius_threshold, this->points); delete this->main_mesh; this->main_meshed = false; } template Mesh & Domain::getAuxiliaryMesh() { if (this->auxiliary_mesh == NULL) { this->initialiseAuxiliaryMesh(); } return *this->auxiliary_mesh; } template Tree & Domain::getTree() { return *this->tree; } template bool Domain::inFemDomain(const PointRef & point) { return this->main_mesh->containsExactNode(point); } template void Domain::fill_atoms() { this->atoms.clear(); for (Uint point_id = 0; point_id < this->points.getSize() ; ++point_id) { PointRef point = this->points.getPoint(point_id); if (this->refined->is_inside(point, this->tol)) { this->atoms.addPoint(point); } } } template void Domain::fill_coupling_atoms() { this->coupling_atom_ids.clear(); for (Uint atom_id = 0 ; atom_id < this->atoms.getSize() ; ++atom_id) { const PointRef & point = this->atoms.getConstPoint(atom_id); if (!this->atomistic_skinned->is_inside(point, this->tol) && this->inFemDomain(point)) { this->coupling_atom_ids.push_front(atom_id); } } } template void Domain::dump_atoms_lammps(std::string filename) { if (this->atoms.getSize() == 0) { throw(std::runtime_error("there are no atoms")); } this->atoms.dump_lammps(filename); } template void Domain::dump_atoms_paraview(std::string filename) { if (this->atoms.getSize() == 0) { throw(std::runtime_error("there are no atoms")); } this->atoms.dump_paraview(filename); } template const PointContainer & Domain::getAtoms() { if (this->atoms.getSize() == 0) { throw(std::runtime_error("there are no atoms. Did you forget to call Domain::fill_atoms() on this domain?")); } return this->atoms; } template bool Domain::inDomain(PointRef & point) { InsideObject eval = this->overall->from_border(point, this->tol); if (!eval.is_inside) { // we're out of the domain return false; } else { if (this->boundaries_special_treatment) { int int_refinement = this->interpolate(point); Real refinement = int_refinement * this->lattice->getRepresentativeConstant(); return max(-eval.distance, 0.) >= refinement*.75; // heuristic value } else { return true; } } } template void Domain::compute_coupling_atoms(PointContainer & interface, PointContainer & pad, const Geometry & atomistic) { if (this->atoms.getSize() == 0) { throw("there are no atoms"); } + if (!this->main_meshed) { + throw std::runtime_error("this domain has not been meshed. Did you call getMesh(quality)?"); + } if (&atomistic == NULL) { - this->getMesh().compute_interface_atoms(interface, *this->atomistic, + this->main_mesh->compute_interface_atoms(interface, *this->atomistic, this->tol); } else { - this->getMesh().compute_interface_atoms(interface, atomistic, + this->main_mesh->compute_interface_atoms(interface, atomistic, this->tol); } // now, clean the pad atoms to get rid of the interface atoms for (Uint interface_index = 0 ; interface_index < interface.getSize() ; ++interface_index) { PointRef interface_point = interface.getPoint(interface_index); auto pad_index = this->coupling_atom_ids.begin(); auto pad_before_index = this->coupling_atom_ids.before_begin(); auto pad_end = this->coupling_atom_ids.end(); bool searching = true; for (; pad_index != pad_end && searching; ++pad_index, ++pad_before_index) { PointRef pad_point = this->atoms.getPoint(*pad_index); if (pad_point == interface_point) { this->coupling_atom_ids.erase_after(pad_before_index); searching=false; } } } auto pad_index = this->coupling_atom_ids.begin(); auto pad_end = this->coupling_atom_ids.end(); for(;pad_index != pad_end; ++pad_index) { pad.addPoint(this->atoms.getPoint(*pad_index)); } } +/* -------------------------------------------------------------------------- */ +template +void Domain::insertCustomPoint(const PointRef & point) { + this->points.addPoint(point); + if (this->main_meshed) { + throw std::runtime_error("This domain was aleady meshed. The new point was inserted, but may just have fucked up the domain" + ); + } +} + template class Domain; template class Domain; diff --git a/src/domain/domain.hh b/src/domain/domain.hh index d423f56..0bf67d6 100644 --- a/src/domain/domain.hh +++ b/src/domain/domain.hh @@ -1,236 +1,239 @@ /** * @file domain.hh * @author Till Junge * @date Thu Apr 19 17:34:16 2012 * * @brief defines the domain class. * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #ifndef __CADD_MESHER_DOMAIN_HH__ #define __CADD_MESHER_DOMAIN_HH__ #include "common.hh" #include "cadd_mesh.hh" #include "lattice.hh" #include "point_container.hh" #include #include #include template class Tree; template class Domain { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /*! \param overall_ computational domain * \param refined_ subset of \a overall_ where atomistic resolution is required * in the case of CADD, this includes pad and MD * \param atomistic_ subset of \a refined_. Contains all points which will be * real atoms, i.e., interface and free * \param atomistic_skinned_ subset of \a atomistic_. contains all free non- interface atoms * \param lattice_ * \param centre The origin from where the domain is to be filled * \param boundaries_special_treatment whether the boundaries should be meshed * seperately. This is necessary if the overall sim box size cannot be a * power-of-2-multiple of a lattice or if the lattice is not aligned with * the overall geometry * \param tolerance when deciding whether a point is inside the domain or not, - * this tol is used inclusively on the one side and exclusively on the other + * this tol is used inclusively on the one side and exclusively on the other: + * positive tol is exclusive on the upper side and inclusive on the lower */ Domain(const Geometry & overall_, const Geometry & refined_, const Geometry & atomistic_, const Geometry & atomistic_skinned_, const Lattice & lattice_, - bool boundaries_special_treatment = 0, + bool boundaries_special_treatment = false, const PointRef & centre=*reinterpret_cast *>(NULL), Real tolerance = 1e-4, const std::string & qh_flags = default_flags ); virtual ~Domain(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /*! adds a point-refinement data point. Make sure that in the end, the convex * hull of the points encompasses the entire Geometry \a overall_. * \param coords * \param refinement */ void setPointRefinement(const Real * coords, Real & refinement); void setPointRefinement(const PointRef & coords, Real & refinement){ this->setPointRefinement(coords.getArray(), refinement); } /// fills the domain void fill_points(bool periodicity[DIM]) throw(std::string); void fill_points() throw(std::string) { bool periodicity[DIM] = {0}; this->fill_points(periodicity); } void fill_points(bool period_x, bool period_y, bool period_z = false) throw(std::string) { if ((DIM==twoD) && period_z) { throw("you cannot specify a periodicity in z for a 2D problem"); } bool period[DIM] = {period_x, period_y}; if (DIM == threeD) { period[DIM-1] = period_z; } return this->fill_points(period); } void complement_periodically(int direction, Uint dim) throw(std::string); /* returns whether a point is to be considered for the fem mesh * \param point*/ bool inFemDomain(const PointRef & point); bool inDomain(PointRef & point); /// simple exposure of Geometry::from_border inline Real from_border(const Real * point, const Real & tol) const { return this->overall->from_border(point, tol).distance;} /*! returns the interpolated refinement at point \a coords * \param coords */ inline Real interpolate(const Real * coords) const { return this->auxiliary_mesh->interpolate(coords, this->refinement);} inline Real interpolate(PointRef point) const { return this->interpolate(&(point.getComponent(0)));} /*! dumps the atoms in lammps format * \param filename */ void dump_atoms_lammps(std::string filename); /*! dumps the atoms in paraview format * \param filename */ void dump_atoms_paraview(std::string filename); /*! computes the pad and interface atoms. * Theoretically, the atomistic geometry does not need to be specified. In * Practice however, it seems to be important frequently to fiddle around * with it in order to deal correctly with periodic boundary conditions. */ void compute_coupling_atoms(PointContainer & interface, PointContainer & pad, const Geometry & atomistic = *reinterpret_cast* >(NULL)); /*! Deprecated legacy function. Use compute_coupling_atoms instead */ void compute_interface_atoms(PointContainer & interface, PointContainer & pad, const Geometry & atomistic = *reinterpret_cast* >(NULL)) { std::cout << "Deprecation warning. Use compute_coupling_atoms instead of " << "compute_interface_atoms!" << std::endl; this->compute_coupling_atoms(interface, pad, atomistic);} + /// use with care + void insertCustomPoint(const PointRef & point); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// returns a pointer to the lattice inline Lattice * getLattice() {return this->lattice;} /// returns a reference to the points inline const PointContainer & getPoints() const {return this->points;} /*! returns a reference to the main mesh. If the mesh does not yet exist, it * is created and then returned */ Mesh & getMesh(Real quality = 10) throw(std::string); /*! returns a reference to the auxiliary mesh. Throws an error \a std::string * if the mesh does not yet exist */ Mesh & getAuxiliaryMesh(); Tree & getTree(); const PointContainer & getAtoms(); void splitDegenerates(const Real & radius_threshold); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// Creates the auxiliary Mesh using the points defined in setPointRefinement void initialiseAuxiliaryMesh(); void fill_atoms(); void fill_coupling_atoms(); // geometric definitions /// computational domain Geometry * overall; /// subdomain of \a overall, atomistic refinement here Geometry * refined; /// subdomain of \a refined. all real atoms, no pad Geometry * atomistic; /// subdomain of \a atomistic. real atoms without interface atoms Geometry * atomistic_skinned; // lattice Lattice * lattice; // refinement definition /// points used as nodes for the auxiliary mesh PointContainer refinement_points; /// nodal refinement values to be interpolated std::vector refinement; /// auxiliary mesh, used to compute refinement in any point of the computational domain Mesh * auxiliary_mesh; /// quadtree (2D) or octtree (3D) used to make a continuously refining mesh Tree * tree; // points definition /// all points, including pure atoms and pure fem nodes PointContainer points; //this is the main container with nodes and atoms /// atoms, including pad PointContainer atoms; /// pad and interface atom indices. this is filled when the atoms are /// computed, so that we can separate interface and pad std::forward_list coupling_atom_ids; /// nodes, including interface atoms PointContainer nodes; /// main mesh object Mesh * main_mesh; /// whether \a points has been filled bool filled; /// whether \a main_mesh has been meshed bool main_meshed; /*! whether the boundaries should be meshed seperately (important if the * overall size cannot be a multiple of 2 lattices in every direction. */ bool boundaries_special_treatment; Real mesh_quality; /*! tolerance to be used to decide whether a point is within the domain This is an inclusive tolerance on the lower bounds of the domain and exlusive at the upper bounds */ Real tol; std::string qh_flags; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator <<(std::ostream & stream, const Domain & _this) { _this.printself(stream); return stream; } #endif /* __CADD_MESHER_DOMAIN_HH__ */ diff --git a/src/domain/mini_mesh.cc b/src/domain/mini_mesh.cc index e7cddff..314de88 100644 --- a/src/domain/mini_mesh.cc +++ b/src/domain/mini_mesh.cc @@ -1,118 +1,144 @@ #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; + if (std::isnan(this->grad.at(DIM*local_nd_id+i))) { + std::stringstream error; + error << "caught a NaN at grad computation" << std::endl + <<"(new_quality-ref_quality) = " << (new_quality-ref_quality) + << std::endl << "step size = " << step_size << std::endl + << "(new_quality-ref_quality)/step_size" + << (new_quality-ref_quality)/step_size << std::endl; + throw NanStepError(error.str()); + } } } } /* -------------------------------------------------------------------------- */ 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])); } + if (max_step == 0) { + std::stringstream error("gradient is zero!"); + throw ZeroStepError(error.str()); + } Real factor = step_size/max_step; for (Real & val: this->grad) { val *= factor; + if (std::isnan(val)) { + std::stringstream error; + error << "caught a NaN after normalisation. Normalisation factor is " + << factor; + this->updateGradient(step_size); + this->normGradient(step_size); + throw NanStepError(error.str()); + } + } +} + +/* -------------------------------------------------------------------------- */ +template +void MiniMesh::stepGradient(Real factor) { + Uint nb_free_nodes = this->free_nodes.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] += factor*this->grad[DIM*local_nd_id+i]; + } } } /* -------------------------------------------------------------------------- */ 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]; - } + try { + this->normGradient(initial_step_size); + this->stepGradient(-1.); + } catch (const ZeroStepError & e) { } + old_quality = quality; quality = this->computeQuality(); // make sure it was an improvement - if ((old_quality < quality) && (initial_step_size > min_step_size)) { + 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]; - } - } + this->stepGradient(1.); 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 index 2bd97ed..95af577 100644 --- a/src/domain/mini_mesh.hh +++ b/src/domain/mini_mesh.hh @@ -1,49 +1,72 @@ #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(){}; }; + class FiniteStepError: public std::runtime_error { + public: + FiniteStepError(std::string str):std::runtime_error(str){}; + virtual ~FiniteStepError(){}; + }; + + class NanStepError: public FiniteStepError { + public: + NanStepError(std::string str):FiniteStepError(str){}; + virtual ~NanStepError(){}; + }; + + class ZeroStepError: public FiniteStepError { + public: + ZeroStepError(std::string str):FiniteStepError(str){}; + virtual ~ZeroStepError(){}; + }; + + 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); + /// step along the gradient + void stepGradient(Real factor = 1.); + 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 262fca3..6b298e2 100644 --- a/src/geometry/cylinder.cc +++ b/src/geometry/cylinder.cc @@ -1,388 +1,387 @@ #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); + keep_max(violation, -dot_p-tol); + keep_max(violation, dot_p+tol - 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] = { 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] = {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] = {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(); + tol = .1;//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;