diff --git a/src/domain/cadd_mesh.cc b/src/domain/cadd_mesh.cc index a77d057..29daf80 100644 --- a/src/domain/cadd_mesh.cc +++ b/src/domain/cadd_mesh.cc @@ -1,916 +1,916 @@ 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" 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-3); + (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 * 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); } } 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); 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)); 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 class Mesh<2>; template class Mesh<3>; diff --git a/src/lattice/lattice.cc b/src/lattice/lattice.cc index 4af12b8..ce97b02 100644 --- a/src/lattice/lattice.cc +++ b/src/lattice/lattice.cc @@ -1,239 +1,260 @@ /** * @file lattice.cc * @author Till Junge * @date Fri Apr 13 15:55:56 2012 * * @brief * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #include "lattice.hh" #include #include "gsl/gsl_vector.h" #include "gsl/gsl_matrix.h" #include #include #include #include template Real determinant(gsl_matrix * rotation) throw(std::string){ Real determinant = 0; Uint num_terms; if (DIM == 2) { num_terms = 1; } else if (DIM == 3) { num_terms = 3; } else { std::stringstream error_stream ; error_stream << "determinant for matrices other than 2x2 or 3x3 are not " << "implemented" << std::endl; std::cout << error_stream.str(); throw(error_stream.str()); } for (Uint i = 0; i < num_terms ; ++i) { Real increment = 1; Real decrement = 1; for (Uint j = 0; j < DIM ; ++j) { increment *= gsl_matrix_get(rotation, j, (j+i)%DIM); decrement *= gsl_matrix_get(rotation, DIM-1-j, (j+i)%DIM); } determinant += increment; determinant -= decrement; } return determinant; } template Lattice::Lattice(Real * _constants, Real * miller) throw(std::string) :atoms("atoms") { this->constants[DIM] = 1; for (Uint i = 0; i < DIM; ++i) { this->constants[i] = _constants[i]; this->constants[DIM] *= _constants[i]; } this->constants[DIM] = pow(this->constants[DIM], 1/static_cast(DIM)); this->millers = gsl_matrix_calloc(DIM, DIM); this->rotation_matrix = gsl_matrix_calloc(DIM, DIM); if (miller == NULL) { gsl_vector_view diag = gsl_matrix_diagonal(this->millers); gsl_vector_set_all(&diag.vector, 1); diag = gsl_matrix_diagonal(this->rotation_matrix); gsl_vector_set_all(&diag.vector, 1); } else { for (Uint i = 0; i < DIM; ++i) { for (Uint j = 0; j < DIM; ++j) { gsl_matrix_set(this->millers, i, j, miller[i*DIM + j]); } gsl_vector_const_view miller_row = gsl_matrix_const_row(this->millers, i); Real norm = gsl_blas_dnrm2(&miller_row.vector); gsl_vector_view rotation_column = gsl_matrix_column(this->rotation_matrix, i); gsl_vector_memcpy(&rotation_column.vector, &miller_row.vector); gsl_vector_scale(&rotation_column.vector, 1./norm); } } Real det = determinant(this->rotation_matrix); if (fabs(det - 1) > std::numeric_limits::epsilon()){ std::stringstream error_stream; error_stream << "Error: the miller indices you specified yield an invalid " << "rotation matrix (The determinant is " << det << " instead " << "of 1)." << std::endl; for (Uint i = 0; i < DIM ; ++i) { error_stream << "["; for (Uint j = 0; j < DIM-1 ; ++j) { error_stream << gsl_matrix_get(this->rotation_matrix, i, j) << ", "; } error_stream << gsl_matrix_get(this->rotation_matrix, i, DIM-1) << "]" << std::endl; } error_stream <<" Make sure they for a normal basis and comply with " << "the right hand rule. Lattice details: " << std::endl; this->printself(error_stream, indent_incr); std::cout << error_stream.str(); throw error_stream.str(); } // compute spatial constants this->spatial_constants[DIM] = 1; for (Uint dir = 0 ; dir < DIM ; ++dir) { this->spatial_constants[dir] = 0; for (Uint j = 0 ; j < DIM ; ++j) { this->spatial_constants[dir] += square(this->constants[j] * gsl_matrix_get(this->millers, dir, j)); } this->spatial_constants[dir] = sqrt(spatial_constants[dir]); this->spatial_constants[DIM] *= this->spatial_constants[dir]; } this->spatial_constants[DIM] = pow(this->spatial_constants[DIM], 1/static_cast(DIM)); } template Lattice::Lattice(const Lattice & other) { for (Uint i = 0 ; i < DIM+1 ; ++i) { this->constants[i] = other.constants[i]; this->spatial_constants[i] = other.spatial_constants[i]; } this->millers = gsl_matrix_calloc(DIM, DIM); gsl_matrix_memcpy (this->millers, other.millers); this->rotation_matrix = gsl_matrix_calloc(DIM, DIM); gsl_matrix_memcpy(this->rotation_matrix, other.rotation_matrix); this->atoms = other.atoms; } template Lattice::~Lattice() { gsl_matrix_free(this->millers); gsl_matrix_free(this->rotation_matrix); } + + +/* -------------------------------------------------------------------------- */ template void Lattice::fillAtoms(int * lattice_coords, PointContainer & container) { gsl_vector * current_point = gsl_vector_alloc(DIM); gsl_vector * rotated_point = gsl_vector_alloc(DIM); for (Uint i = 0; i < this->atoms.getSize(); ++i) { PointRef current_atom = this->atoms.getPoint(i); for (Uint j = 0; j < DIM; ++j) { gsl_vector_set(current_point, j, current_atom.getComponent(j) + lattice_coords[j] * this->constants[j]); } gsl_blas_dgemv (CblasNoTrans, 1.0, this->rotation_matrix, current_point, 0.0, rotated_point); container.addPoint(rotated_point->data); } gsl_vector_free(current_point); gsl_vector_free(rotated_point); } +/* -------------------------------------------------------------------------- */ +template +void Lattice::fillAtoms(const std::vector sw_coords, + PointContainer & container) { + if (sw_coords.size() != DIM) { + std::stringstream error; + error << "expected a vector of dimension " << DIM << " but got " + << sw_coords.size() << "."; + std::runtime_error(error.str()); + } + Real coords[DIM]; + for (Uint i = 0 ; i < DIM ; ++i) { + coords[i] = sw_coords[i]; + } + this->fillAtoms(coords, container); +} + +/* -------------------------------------------------------------------------- */ template void Lattice::fillAtoms(Real * sw_coords, PointContainer & container) { gsl_vector * current_point = gsl_vector_alloc(DIM); gsl_vector * rotated_point = gsl_vector_alloc(DIM); for (Uint i = 0; i < this->atoms.getSize(); ++i) { PointRef current_atom = this->atoms.getPoint(i); for (Uint j = 0; j < DIM; ++j) { gsl_vector_set(current_point, j, current_atom.getComponent(j)); } gsl_blas_dgemv (CblasTrans, 1.0, this->rotation_matrix, current_point, 0.0, rotated_point); for (Uint i = 0; i < DIM ; ++i) { rotated_point->data[i] += sw_coords[i]; } container.addPoint(rotated_point->data); } gsl_vector_free(current_point); gsl_vector_free(rotated_point); } template const Real & Lattice::getMinConstant() const { Real tmp = std::numeric_limits::max(); Uint index = 0; for (Uint dir = 0 ; dir < DIM ; ++dir) { if (tmp < this->constants[dir]) { tmp = this->constants[dir]; index = dir; } } return this->constants[index]; } template void Lattice::printself(std::ostream & stream, int indent) const { indent_space(stream, indent); stream << this->lattice_type << DIM << "D Lattice, a = ("; for (Uint i = 0; i < DIM-1; ++i) { stream << this->constants[i] << ", "; } stream << this->constants[DIM-1] << "):" << std::endl; this->atoms.printself(stream, indent+indent_incr); indent_space(stream, indent+indent_incr); stream << "Miller orientations :" << std::endl; for (Uint i = 0; i < DIM; ++i) { indent_space(stream, indent + 2*indent_incr); stream << "["; for (Uint j = 0; j < DIM-1; ++j) { stream << gsl_matrix_get(this->millers, i, j) << ", "; } stream << gsl_matrix_get(this->millers, i, DIM-1) << "]" << std::endl; } stream << "Transposed Rotation matrix :" << std::endl; for (Uint i = 0; i < DIM; ++i) { indent_space(stream, indent + 2*indent_incr); stream << "["; for (Uint j = 0; j < DIM-1; ++j) { stream << gsl_matrix_get(this->rotation_matrix, i, j) << ", "; } stream << gsl_matrix_get(this->rotation_matrix, i, DIM-1) << "]" << std::endl; } } template std::string Lattice::lattice_type = ""; template class Lattice<2>; template class Lattice<3>; diff --git a/src/lattice/lattice.hh b/src/lattice/lattice.hh index e367c0f..c7673d4 100644 --- a/src/lattice/lattice.hh +++ b/src/lattice/lattice.hh @@ -1,134 +1,137 @@ /** * @file lattice.hh * @author Till Junge * @date Fri Apr 13 11:08:40 2012 * * @brief interface to lattices * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #ifndef __CADD_MESHER_LATTICE_HH__ #define __CADD_MESHER_LATTICE_HH__ #include "common.hh" #include "point_container.hh" #include "gsl/gsl_matrix.h" #include +#include /*! \class Lattice Base class of lattices, can also be used to create custom lattices */ template class Lattice { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /*! Constructor * \param constants array of lattice constants in \a x, \a y and \a z * direction of the lattice * \param miller matrix of Miller indices (as 1D array of size \a DIM * \a DIM * ). The rows are the Miller vectors to be aligned with the corresponding * direction of the global coordinate system. If \a NULL, the identity is * used.*/ Lattice(Real * constants, Real * miller = NULL) throw(std::string); Lattice(const Lattice & other); virtual ~Lattice(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual Lattice * resolveType() const = 0; /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; void fillAtoms(int * lattice_coords, PointContainer & container); /*! fill the atoms into a PointContainer with an offset given * by the coordinates fo the south-west corner of the lattice * \param sw_coords coordinates of the south-west corner of the lattice * \param container reference to PointContainer to fill */ + void fillAtoms(const std::vector sw_coords, + PointContainer & container); void fillAtoms(Real * sw_coords, PointContainer & container); /*! returns the rotation matrix which rotates the lattice back into reference * configuration*/ gsl_matrix * getRotationTranspose() { return this->rotation_matrix;} /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /*! returns the lattice constant in a given direction (in the lattice * reference orientation) * \param direction */ inline const Real & getConstant(Uint direction) const {return this->constants[direction];} /*! returns the geometric mean of the constants. Hopefully this is a good measure. * The geometric mean is defined as \f[G = \sqrt[n]{x_1 x_2 \cdots x_n}\f] */ inline const Real & getRepresentativeConstant() const {return this->constants[DIM];} /*! computes the size of a periodic cell of the rotated lattice in the global coordinate sytem i.e: - a lattice with Miller indices + a lattice with Miller indices [ -1, 1, 0] M = [ 1, 1, 1] , [ 1, 1, -2] where the rows are the Miller vectors, the spatial constants (for a cubic lattice) are sqrt(2), sqrt(3) and sqrt(6) in x, y, z respectively */ inline const Real & getSpatialConstant(Uint direction) const { return this->spatial_constants[direction];} /*! returns the array of lattice constants */ inline Real * getConstants() {return this->constants;} inline Real * getSpatialConstants() {return this->spatial_constants;} const Real & getMinConstant() const; /*! returns the direction of a lattice vector in the global coordinate system * this is a unit vector * \param direction */ inline const Real * getLatticeVector(Uint direction) { return gsl_matrix_ptr(this->rotation_matrix, direction, 0); } inline std::vector getLatticeStdVector(Uint direction) { const Real * latt_vec = this->getLatticeVector(direction); return std::vector(latt_vec, latt_vec + DIM); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: Real constants [DIM + 1]; gsl_matrix * millers; gsl_matrix * rotation_matrix; PointContainer atoms; Real spatial_constants[DIM + 1]; private: static std::string lattice_type; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator <<(std::ostream & stream, const Lattice & _this) { _this.printself(stream); return stream; } template Lattice * newLattice(const subLattice & lattice) { return static_cast*>(new subLattice(lattice)); } #endif /* __CADD_MESHER_LATTICE_HH__ */ diff --git a/test/py_wrapper_test.py.cmake b/test/py_wrapper_test.py.cmake index 6a33121..ce4b5aa 100644 --- a/test/py_wrapper_test.py.cmake +++ b/test/py_wrapper_test.py.cmake @@ -1,340 +1,332 @@ #!/usr/bin/env python from __future__ import print_function from __future__ import division import sys sys.path.append('@CADD_MESH_PY_BUILD_PATH@') import math import numpy as np from CADDMesher import * error_counter = 0 def vector_test(): a = vecReal([1,2,3,4,5]) #print("vector a = [{0}]".format(", ".join([str(item) for item in a]))) return False def fcc_lattice_test(): consts = vecReal([1, 3, 4]) miller = vecReal([1, 0, 0, 0, 0, 1, 0, 1, 0]) lattice = FccLattice2d(consts, miller) rep_const = lattice.getRepresentativeConstant() ref_val = 3**(1/2) if rep_const != ref_val: raise Exception( "Representative Constant wrong: it's {0} but should be {1}\n".format( rep_const, ref_val) + "difference = {0}".format(ref_val - rep_const)) return False def block_test(): b = Block2d([0,1,0,2], 'chicken') return False def PointRef_test(): b = PointRef2d(vecReal([1,2])) c = PointRef3d(vecReal([1,2,3])) d = PointRef2d(b) return False def PointContainer_test(): b = PointRef2d(vecReal([1,2])) cont = PointContainer2d("MonkeyContainer") cont.addPoint(b) cont.addPoint(vecReal([12.4, 15.2, 83.1])) point = cont.getPoint(0) return False def Mesh_test(): lim = 32 points = list() nodal_vals = list() for i in range(lim+1): for j in range(lim+1): points.append(i/lim) points.append(j/lim) nodal_vals.append((i+j)/lim) minmax = [.45, 10.6, .45, .55] exclude = Block2d(minmax, "excludes some stuff") minmax = [.4, 10.6, .4, .6] atomic = Block2d(minmax, "atomistic domain") mesh = Mesh2d(points, 0, exclude) mesh.delaunay(1000) mesh.cleanup(atomic) mesh.dump_paraview("after_cleanup_2D_python") return False def Domain_test(): consts = vecReal([4.04, 4.04]) miller = vecReal([ 1., 1., -1., 1.]) lattice = FccLattice2d(consts, miller) ## print("Lattice = {0}".format(lattice)) # Prepare the domain geometries min_max = vecReal([-200, 200, -300, 0]) overall = Block2d(min_max, "overall"); min_max = vecReal([-25, 25, -35, 0]) refined = Block2d(min_max, "refined") min_max = vecReal([-20, 20, -30, 0]) atomistic = Block2d(min_max, "atomistic") min_max = vecReal([-19, 19, -29, 0]) atomistic_skinned = Block2d(min_max, "atomistic_skinned") # Create domain domain = Domain2d(overall, refined, atomistic, atomistic_skinned, lattice, True) refinement_point = PointRef2d() punkte = PointContainer2d("just for debugging") # set up auxiliary mesh for i in range(2): for k in range(2): refinement_point[0] = (1-2*i)*201 refinement_point[1] = -(k * 301 - (1-k)) refinement = float(8); domain.setPointRefinement(refinement_point, refinement) punkte.addPoint(refinement_point) refinement_point[0] = (1-2*i)*25 refinement_point[1] = -(k * 35 - (1-k)) refinement = 0.; domain.setPointRefinement(refinement_point, refinement) punkte.addPoint(refinement_point) refinement_point[0] = (1-2*i)*100 refinement_point[1] = -(k * 130 - (1-k)) refinement = 4.; domain.setPointRefinement(refinement_point, refinement) punkte.addPoint(refinement_point) ## print(punkte) testmesh = Mesh2d(punkte.getVector(), 0) ## print("Test mesh before delaunay = \n{0}".format(testmesh)) testmesh.delaunay(10000); testmesh.initInterpolation() ## print("Test mesh = \n{0}".format(testmesh)) domain.fill_points(); print("Number of tree boxes = {0}".format(domain.getTree().getNbBoxes())) domain.getMesh().dump_paraview("domain_test_output_2d_python") print("domain.getPoint() : \n{0}".format(domain.getPoints())) return False def Domain_test_3d(): consts = vecReal([4.04, 4.04, 4.04]) miller = vecReal([ 1, 1, 0, -1, 1, 0, 0, 0, 1]) lattice = FccLattice3d(consts, miller) # Prepare the domain geometries min_max = vecReal([-200, 200, -200, 200, -300, 0]) overall = Block3d(min_max, "overall"); min_max = vecReal([-25, 25, -25, 25, -35, 0]) refined = Block3d(min_max, "refined") min_max = vecReal([-20, 20, -20, 20, -30, 0]) atomistic = Block3d(min_max, "atomistic") min_max = vecReal([-19.9, 19.9, -19.9, 19.9, -29.9, 0]) atomistic_skinned = Block3d(min_max, "atomistic_skinned") # Create domain domain = Domain3d(overall, refined, atomistic, atomistic_skinned, lattice, True) refinement_point = PointRef3d() # set up auxiliary mesh for i in range(2): for j in range(2): for k in range(2): refinement_point[0] = (1-2*i)*201 refinement_point[1] = (1-2*j)*201 refinement_point[2] = -(k * 301 - (1-k)) print(refinement_point) refinement = float(8); domain.setPointRefinement(refinement_point, refinement) refinement_point[0] = (1-2*i)*25 refinement_point[1] = (1-2*j)*25 refinement_point[2] = -(k * 35 - (1-k)) refinement = 0.; domain.setPointRefinement(refinement_point, refinement) refinement_point[0] = (1-2*i)*100 refinement_point[1] = (1-2*j)*100 refinement_point[2] = -(k * 130 - (1-k)) refinement = 4.; domain.setPointRefinement(refinement_point, refinement) domain.fill_points(); print("Number of tree boxes = {0}".format(domain.getTree().getNbBoxes())) return False def Cylinder_test_3d(): consts = vecReal([4.04, 4.04, 4.04]) miller = vecReal([ 1, 1, 0, -1, 1, 0, 0, 0, 1]) lattice = FccLattice3d(consts, miller) # Prepare the domain geometries base = vecReal([0, 0, 0]) length = 40. axis = vecReal([0, 0, length]) outer_radius = 100 inner_radius = 20 radius = outer_radius overall = Cylinder3d(base, radius, axis, "overall"); radius = inner_radius refined = Cylinder3d(base, radius, axis, "refined") radius = 18 atomistic = Cylinder3d(base, radius, axis, "atomistic") radius = 17 atomistic_skinned = Cylinder3d(base, radius, axis, "atomistic_skinned") # Create domain domain = Domain3d(overall, refined, atomistic, atomistic_skinned, lattice, True) refinement_point = PointRef3d() rx = np.array([1., 0., 0.]) ry = np.array([0., 1., 0.]) # set up auxiliary mesh - - for i in range(10): - angle = 2*math.pi/180*i/10 + angular_resolution = 48 + for i in range(angular_resolution): + angle = 2*np.pi*i/angular_resolution r = np.cos(angle)*rx + np.sin(angle)*ry for k, z in zip([0, -1.], [0, length+1]) : refinement_point[0] = 1.1*outer_radius*r[0] refinement_point[1] = 1.1*outer_radius*r[1] refinement_point[2] = z refinement = float(8); domain.setPointRefinement(refinement_point, refinement) refinement_point[0] = 1.1*inner_radius*r[0] refinement_point[1] = 1.1*inner_radius*r[1] refinement_point[2] = z refinement = float(0); domain.setPointRefinement(refinement_point, refinement) - - domain.fill_points(); - print("Number of tree boxes = {0}".format(domain.getTree().getNbBoxes())) + domain.getAuxiliaryMesh().dump_paraview("auxiliary") + domain.fill_points(); + domain.getMesh().dump_paraview("cylindertest") return False def rotationTest(): a = 4.04 consts = vecReal([a, a, a]) miller = vecReal([-1, 1, 0, 1, 1, 1, 1, 1,-2]) miller = vecReal([ 1, 1, 0, -1, 1, 0, 0, 0, 1]) lattice = FccLattice3d(consts, miller) x_min = y_min = z_min = 0 x_max = math.sqrt(miller[0]**2 + miller[1]**2 + miller[2]**2)*a y_max = math.sqrt(miller[3]**2 + miller[4]**2 + miller[5]**2)*a z_max = math.sqrt(miller[6]**2 + miller[7]**2 + miller[8]**2)*a coord_max = max(abs(x_min), abs(x_max), abs(y_min), abs(y_max), abs(z_min), abs(z_max))*1.1 #Creating a single lattice min_max = vecReal([ x_min, x_max, y_min, y_max, z_min, z_max]) min_max = vecReal([-x_max, x_max, -y_max, y_max, -z_max, z_max]) # Prepare the domain geometries overall = Block3d(min_max, "overall"); # Create domain boundary_special_treatment = False domain = Domain3d(overall, overall, overall, overall, lattice, boundary_special_treatment) refinement_point = PointRef3d() temp_cont = PointContainer3d("temporary container") refinement = 0 # set up auxiliary mesh for i in [-1, 1]: refinement_point[0] = i*coord_max for j in [-1, 1]: refinement_point[1] = j*coord_max for k in [-1, 1]: refinement_point[2] = k*coord_max print(refinement_point) domain.setPointRefinement(refinement_point, refinement) temp_cont.addPoint(refinement_point) domain.fill_points() - domain.getMesh() print("Number of tree boxes = {0}".format(domain.getTree().getNbBoxes())) - print("ro-fucking-tation") - - print("bla") - atoms = PointContainer3d(domain.getAtoms()) - atoms.dump_paraview("rotation_test") - - print("ro-fucking-tation") return False def tester(function, *args, **kwargs): try: ret_val = function(*args, **kwargs) print("Test '{0}' passed".format(function.__name__)) if ret_val != False: raise Exception("Test '{0}' gave return code '{1}'".format( function.__name__, ret_val)) return False except Exception as err: global error_counter error_counter += 1 indented_message = '\n'.join( [' "{0}"'.format(line) for line in err.message.split('\n')]) print('------', err, '--------') print("!!!!! error number {0}".format(error_counter)) print("Test '{0}' failed with the following error:\n{1}".format( function.__name__, indented_message)) return True def main(): overall_failure = False overall_failure |= tester(vector_test) overall_failure |= tester(fcc_lattice_test) overall_failure |= tester(block_test) overall_failure |= tester(PointRef_test) overall_failure |= tester(PointContainer_test) overall_failure |= tester(Mesh_test) overall_failure |= tester(Domain_test) overall_failure |= tester(rotationTest) - #overall_failure |= tester(Cylinder_test_3d) + overall_failure |= tester(Cylinder_test_3d) if not overall_failure: print("\nOverall sucess. All tests passed.") else: print("\nOverall failure, some tests failed!") main()