diff --git a/src/common/common.hh b/src/common/common.hh index a52f65d..f990169 100644 --- a/src/common/common.hh +++ b/src/common/common.hh @@ -1,93 +1,93 @@ /** * @file common.hh * @author Till Junge * @date Wed Apr 4 15:59:37 2012 * * @brief common parts for the cadd_mesher * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #ifndef __cadd_mesh_COMMON_HH__ #define __cadd_mesh_COMMON_HH__ #include #include #include typedef double Real; typedef unsigned int Uint; const static Uint twoD = 2; const static Uint threeD = 3; const static Uint indent_incr = 2; void indent_space(std::ostream & stream, int indent = 0); template inline void keep_max(numeric & value, const numeric & comparison) { if (comparison > value) { value = comparison; } } template inline void keep_min(numeric & value, const numeric & comparison) { if (comparison < value) { value = comparison; } } template void print_vector(const std::vector & vec, const std::string & name, Uint stride = 1, int indent = 0, std::ostream & stream = std::cout){ indent_space(stream, indent); stream << "Vector '" << name << "': " < max_size) { stream.width(7); indent_space(stream, indent+indent_incr); stream << i << " : "; for (Uint j = 0; j < stride-1; ++j) { stream.width(4); stream << std::showpos << "..." << ",\t"; } stream.width(4); stream << "..." << std::endl; } } template -void print_array(el_typo * array, Uint length, const std::string & name) { - std::cout << "Array '" << name << "': (" ; +void print_array(el_typo * array, Uint length, const std::string & name, std::ostream & stream=std::cout) { + stream << "Array '" << name << "': (" ; for (Uint i = 0; i < length-1 ; ++i) { - std::cout << array[i] << ", "; + stream << array[i] << ", "; } - std::cout << array[length-1] << ")" << std::endl; + stream << array[length-1] << ")" << std::endl; } #endif /* __cadd_mesh_COMMON_HH__ */ diff --git a/src/domain/domain.cc b/src/domain/domain.cc index 469f4ae..4c808ad 100644 --- a/src/domain/domain.cc +++ b/src/domain/domain.cc @@ -1,129 +1,136 @@ #include "domain.hh" #include "tree.hh" #include #include template Domain::Domain(Geometry * overall_, Geometry * refined_, Geometry * atomistic_, Geometry * atomistic_skinned_, - Lattice * lattice_, Real * centre) + Lattice * lattice_, Real * centre, + bool boundaries_special_treatment_) :overall(overall_), refined(refined_), atomistic(atomistic_), atomistic_skinned(atomistic_skinned_), lattice(lattice_), auxiliary_mesh(NULL), tree(NULL), main_mesh(NULL), - filled(false), main_meshed(false) + filled(false), main_meshed(false), + boundaries_special_treatment(boundaries_special_treatment_) { - + if (boundaries_special_treatment) { + throw(std::string("The special treatment for the surface is not yet implemented")); + } 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); this->tree = new Tree(centre, 2*pow2size, this->lattice, this); } template Domain::~Domain() { delete this->auxiliary_mesh; delete this->main_mesh; delete this->tree; } template void Domain::setPointRefinement(Real * coords, Real & refinement ) { this->refinement_points.addPoint(coords); this->refinement.push_back(refinement); } template void Domain::fill_points() throw(std::string){ this->initialiseAuxiliaryMesh(); this->tree->fill(this->points); this->filled = true; } 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() { this->auxiliary_mesh = new Mesh(this->refinement_points.getVector()); this->auxiliary_mesh->delaunay(1000); this->auxiliary_mesh->initInterpolation(); } template Mesh & Domain::getMesh() throw(std::string) { if (!this->main_meshed) { this->main_mesh = new Mesh(this->points.getVector(), this->atomistic_skinned); try { this->main_mesh->delaunay(1000); } catch (std::string & error) { this->tree->dump(); this->points.dump(); throw("affen_fehler" + error); } this->main_mesh->cleanup(*this->atomistic); this->main_meshed = true; } return *this->main_mesh; } template Mesh & Domain::getAuxiliaryMesh() throw(std::string) { if (this->auxiliary_mesh == NULL) { throw (std::string("Auxiliary mesh does not exist")); } return *this->auxiliary_mesh; } template Tree & Domain::getTree() { return *this->tree; } template bool Domain::inFemDomain(PointRef & point) { Real distance_from_border = -this->overall->from_border(point); if (distance_from_border < 0) { // we're out of the domain return false; } Real refinement = this->interpolate(point) * this->lattice->getRepresentativeConstant(); bool retval = distance_from_border >= refinement && !this->atomistic_skinned->is_inside(point); return retval; } template bool Domain::inDomain(PointRef & point) { Real distance_from_border = -this->overall->from_border(point); if (distance_from_border < 0) { // we're out of the domain return false; + } else { + if (this->boundaries_special_treatment) { + Real refinement = this->interpolate(point) + * this->lattice->getRepresentativeConstant(); + return distance_from_border >= refinement; + } else { + return true; + } } - Real refinement = this->interpolate(point) - * this->lattice->getRepresentativeConstant(); - - bool retval = distance_from_border >= refinement; - return retval; } template class Domain; template class Domain; diff --git a/src/domain/domain.hh b/src/domain/domain.hh index 93cc323..2162f8a 100644 --- a/src/domain/domain.hh +++ b/src/domain/domain.hh @@ -1,164 +1,173 @@ /** * @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 "mesh.hh" #include "lattice.hh" #include "point_container.hh" #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- * free 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 */ Domain(Geometry * overall_, Geometry * refined_, Geometry * atomistic_, Geometry * atomistic_skinned_, - Lattice * lattice_, Real * centre = NULL); + Lattice * lattice_, Real * centre = NULL, + bool boundaries_special_treatment = false); 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(Real * coords, Real & refinement); /// fills the domain void fill_points() throw(std::string); /* returns whether a point is to be considered for the fem mesh * \param point*/ bool inFemDomain(PointRef & point); bool inDomain(PointRef & point); /// simple exposure of Geometry::from_border inline Real from_border(const Real * point) const { return this->overall->from_border(point);} /*! 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)));} /// returns a pointer to the lattice /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: 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() 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() throw(std::string); Tree & getTree(); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// Creates the auxiliary Mesh using the points defined in setPointRefinement void initialiseAuxiliaryMesh(); // 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; /// 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; + }; /* -------------------------------------------------------------------------- */ /* 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/mesh.cc b/src/domain/mesh.cc index ce2d822..42f57dd 100644 --- a/src/domain/mesh.cc +++ b/src/domain/mesh.cc @@ -1,374 +1,375 @@ extern "C" { #include #include #include #include #include /* for qh_vertexneighbors() */ #include /* for qh_facetarea(facet)*/ } #include "mesh.hh" #include #include #include #include #include #include "dumper_paraview.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); } template Mesh::Mesh(const std::vector & points, Geometry * excluded_points) :permutation(NULL), remutation(NULL), meshed(false), cleaned(false), interpolation_initialised(false) { for (Uint i = 0; i < points.size() ; i += DIM) { if ((excluded_points == NULL) || (!excluded_points->is_inside(&points.at(i))) ) { this->renumerotation.push_back(i/DIM); for (Uint j = 0; j < DIM ; ++j) { this->nodes.push_back(points.at(i+j)); } } } } template Mesh::~Mesh() { if (this->permutation != NULL) { delete this->permutation; delete this->remutation; } } 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, "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); } template void Mesh::delaunay(Real degeneration_criterion){ // make a copy of the nodes before they get jiggled std::vector copy_nodes = this->nodes; boolT ismalloc = False; /* True if qhull should free points in qh_freeqhull() or reallocation */ int exitcode; /* 0 if no error from qhull */ char flags[] = "qhull d QJ Pg Ft"; /* flags for qhull */ /* Call qhull */ exitcode = qh_new_qhull(DIM, copy_nodes.size()/DIM, ©_nodes[0], ismalloc, flags, NULL, stderr); if (exitcode){ std::stringstream exit_stream; exit_stream << "qhull failed with exitcode " << exitcode; throw exit_stream.str(); } else { //I use an stl vector to build the connectivity list. making sure it's empty: connectivity.clear(); //build neighbours: qh_vertexneighbors(); //elements in 3d are 4d facets, hence facetT *facet; vertexT * vertex, ** vertexp; Real element_volume; int nb_elems = 0; std::vector current_element; FORALLfacets { if (facet->upperdelaunay) continue; element_volume = qh_facetarea(facet); ++nb_elems; unsigned int counter = 0; current_element.clear(); FOREACHvertex_(facet->vertices){ current_element.push_back(qh_pointid(vertex->point)); ++ counter; } if (counter != DIM+1){ std::stringstream exit_stream; exit_stream << "Element #" << nb_elems << "has wrong number of vertices: " << counter; throw exit_stream.str(); } else { bool is_valid = degenerate_finder(this->nodes, current_element, element_volume, degeneration_criterion); if (is_valid) { for (unsigned int i = 0; i < DIM+1; ++i) { connectivity.push_back(current_element[i]); } } } } } /* Free the memory */ qh_freeqhull(qh_ALL); this->meshed = true; } template void Mesh::cleanup(const Geometry & atomistic_domain) { Uint nb_points = this->nodes.size()/DIM; 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] += nodes.at(DIM*point_id + k); } if (atomistic_domain.is_inside(&nodes.at(DIM*point_id)) ) { nb_pure_atoms ++; } } // if at least one node is not a pure atom, we're good to go // also, if an element is valid, its nodes are (obviously) valid as well if (nb_pure_atoms < DIM+1) { - for (Uint j = 0; j < DIM+1 ; ++j) { - int point_id = element_nodes[j]; - for (Uint k = 0; k < DIM ; ++k) { - center_gravity[k] /= (DIM+1); - } - if (atomistic_domain.is_inside(center_gravity)) { - break; - } - if (this->permutation->at(point_id) < 0) { - this->permutation->at(point_id) = counter; - this->remutation ->at(counter) = point_id; - ++counter; - for (Uint k = 0; k < DIM ; ++k) { - new_points.push_back(this->nodes.at(point_id*DIM+k)); + // check the center of gravity + for (Uint k = 0; k < DIM ; ++k) { + center_gravity[k] /= (DIM+1); + } + if (!atomistic_domain.is_inside(center_gravity)) { + for (Uint elem_node = 0; elem_node < DIM+1 ; ++elem_node) { + int point_id = element_nodes[elem_node]; + if (this->permutation->at(point_id) < 0) { + this->permutation->at(point_id) = counter; + this->remutation ->at(counter) = point_id; + ++counter; + for (Uint k = 0; k < DIM ; ++k) { + new_points.push_back(this->nodes.at(point_id*DIM+k)); + } } + new_connectivity.push_back(this->permutation->at(point_id)); } - new_connectivity.push_back(this->permutation->at(point_id)); } - } + } } //print_vector(this->nodes, "before cleaning", DIM); this->nodes = new_points; //print_vector(this->nodes, "after cleaning", DIM); //print_vector(*(this->permutation), "permutation"); //print_vector(*(this->remutation), "remutation"); //print_vector(this->connectivity, "before cleaning", DIM+1); this->connectivity = new_connectivity; //print_vector(this->connectivity, "after cleaning", DIM+1); this->cleaned = true; } void inverse(double * A, int N) { int *IPIV = new int[N+1]; int LWORK = N*N; double *WORK = new double[LWORK]; int INFO; dgetrf_(&N,&N,A,&N,IPIV,&INFO); dgetri_(&N,A,&N,IPIV,WORK,&LWORK,&INFO); delete[] IPIV; delete[] WORK; } template 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.at(element_nodes[DIM]*DIM + dir); this->reverse_matrices.push_back(last_node_coord); for (Uint node_local = 0; node_local < DIM ; ++node_local) { this->reverse_matrices.push_back(this->nodes.at(element_nodes[node_local]*DIM + dir) - last_node_coord); } } //invert the matrix to make it useful inverse(&(this->reverse_matrices.at(elem_id*(DIM+1)*(DIM+1))), DIM + 1); } this->interpolation_initialised = true; } template Real Mesh::interpolate(const Real * coords, const std::vector & nodal_values) const { Uint nb_elems = this->connectivity.size()/(DIM+1); Real elem_coords[DIM+2]; //we'll use the last slot for the n-th shape function value Real * shape_vals = elem_coords+1; Real glob_coords[DIM+1]; glob_coords[0] = 1; for (Uint dir = 0; dir < DIM ; ++dir) { glob_coords[dir+1] = coords[dir]; } Real max_viol = -std::numeric_limits::max(); Uint elem_id; bool found = false; int *element_nodes = NULL; //loop through elements and figure out in which one the coords are for (elem_id = 0; elem_id < nb_elems ; ++elem_id) { const Real* matrix = &(this->reverse_matrices.at(elem_id * (DIM+1)*(DIM+1))); cblas_dgemv(CblasRowMajor, CblasNoTrans, DIM+1, DIM+1, 1., matrix, DIM+1, glob_coords, 1, 0., elem_coords, 1); if (fabs(1-elem_coords[0]) > std::numeric_limits::epsilon()){ std::stringstream error_stream ; error_stream << " Reversion matrix for element "<< elem_id << " is obviously wrong, because the sum of " << " Shape functions is not 1 " << std::endl; throw(error_stream.str()); } // check if any of the shape functions are < 0 Real minimum = std::numeric_limits::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 + std::numeric_limits::epsilon() >= 0) { + if (minimum + 2*std::numeric_limits::epsilon() >= 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 = " << std::numeric_limits::epsilon() << std::endl; throw(error_stream.str()); } Real interpolated_value = 0; for (Uint i = 0; i < DIM+1 ; ++i) { int index; if (this->permutation != NULL) { index = this->remutation->at(element_nodes[i]); } else { index = element_nodes[i]; } index = this->renumerotation.at(index); interpolated_value += nodal_values.at(index) * shape_vals[i]; } return interpolated_value; } template void Mesh::dump(std::string filename) throw (std::string){ if (this->meshed) { iohelper::DumperParaview dumper; dumper.SetMode(iohelper::TEXT); dumper.SetPoints(&this->nodes[0], DIM, this->nodes.size()/DIM, filename); iohelper::ElemType type; if (DIM == 2) { type = iohelper::TRIANGLE1; } else if (DIM == 3) { type = iohelper::TETRA1; } else { throw(std::string("Only 2 and 3D meshes can be dumped")); } dumper.SetConnectivity(&this->connectivity[0], type, connectivity.size()/(DIM+1), iohelper::C_MODE); dumper.SetPrefix("./"); dumper.Init(); dumper.Dump(); } else { throw(std::string("this mesh has not been delaunay'ed yet")); } } template class Mesh<2>; template class Mesh<3>; diff --git a/src/domain/tree.cc b/src/domain/tree.cc index 535f76b..ef8c7a4 100644 --- a/src/domain/tree.cc +++ b/src/domain/tree.cc @@ -1,237 +1,236 @@ #include "tree.hh" #include "common.hh" //This function sets a vec1's components to vec2's components template inline void vector_set(Real * vec1, const Real * vec2) { for (Uint i = 0; i < DIM ; ++i) { vec1[i] = vec2[i]; } } // stretches a vector template inline void vector_stretch(Real * vector, const Real * factors) { for (Uint i = 0; i < DIM ; ++i) { vector[i] *= factors[i]; } } //This function adds factor times addition to vec1 template inline void vector_incr(Real * vec1, const Real * addition, const Real & factor) { for (Uint i = 0; i < DIM ; ++i) { vec1[i] += factor * addition[i]; } } template Tree::Tree(Real * centre_, Uint size_, Lattice * lattice_, Domain * domain_) :size(size_), domain(domain_), lattice(lattice_), neighbours(0), level(0) { this->nb_boxes ++; if (centre_ == NULL) { for (Uint i = 0; i < DIM ; ++i) { this->south_west[i] = 0; } } else { vector_set(this->south_west, centre_); } for (Uint latt_dir = 0; latt_dir < DIM ; ++latt_dir) { vector_incr(this->south_west, this->lattice->getLatticeVector(latt_dir), -0.5*this->size * this->lattice->getConstant(latt_dir)); } } template Tree::~Tree() { for (typename std::vector *>::iterator child= this->children.begin(); child != this->children.end(); ++child) { delete *child; } } template void Tree::fill(PointContainer & points) { // the first one always has to be split manually, because none of its corners are // within el domaino if (doIHaveToSplit()) { this->split(); typename std::vector * >::iterator child= this->children.begin(); for (; child != this->children.end() ; ++child) { (*child)->fill(points); } } else { this->create_points(points); } } template void Tree::dump(std::string file_name) { std::ofstream outfile; outfile.open(file_name.c_str()); this->dump(outfile); outfile.close(); } template Tree::Tree(Real * south_west_, Uint size_, Lattice * lattice_, Domain * domain_, Uint neighbours_) :size(size_), domain(domain_), lattice(lattice_), neighbours(neighbours_) { this->nb_boxes ++; for (Uint i = 0; i < DIM ; ++i) { this->south_west[i] = south_west_[i]; } } template void Tree::dump(std::ostream & stream) { Uint nb_corners = 1<(coords, this->south_west); vector_set(scaled_south_west, this->south_west); for (Uint latt_dir = 0; latt_dir < DIM ; ++latt_dir) { Uint factor = reorder[corner_id]>>latt_dir & 1; vector_incr(coords, this->lattice->getLatticeVector(latt_dir), factor * this->size * this->lattice->getConstant(latt_dir)); } for (Uint dir = 0; dir < DIM ; ++dir) { stream << coords[dir] << "\t"; } stream << std::endl; } for (Uint dir = 0; dir < DIM; ++dir) { stream << scaled_south_west[dir] << "\t"; } stream << std::endl << std::endl; typename std::vector * >::iterator child= this->children.begin(); for (; child != this->children.end() ; ++child) { (*child)->dump(stream); } } template void Tree::split() { Uint new_size = size >> 1; Uint geometric_size = new_size; keep_max(geometric_size, static_cast(1)); Uint nb_children = 1< * child = new Tree(this->south_west, new_size, this->lattice, this->domain, this->neighbours); child->level = this->level+1; this->children.push_back(child); } else { for (Uint child_id = 0; child_id < nb_children ; ++child_id) { Real child_sw[DIM]; vector_set(child_sw, this->south_west); for (Uint latt_dir = 0; latt_dir < DIM ; ++latt_dir) { Uint factor = child_id >> latt_dir & 1; vector_incr(child_sw, this->lattice->getLatticeVector(latt_dir), factor * geometric_size*this->lattice->getConstant(latt_dir)); } Tree * child = new Tree(child_sw, new_size, this->lattice, this->domain, ~(this->neighbours | child_id)); this->children.push_back(child); child->level = this->level+1; } } } template bool Tree::doIHaveToSplit() { Uint nb_corners = 1<size == 0 ) { return false; } for (Uint corner_id = 0; corner_id < nb_corners ; ++corner_id) { vector_set(coords, this->south_west); for (Uint latt_dir = 0; latt_dir < DIM ; ++latt_dir) { Uint factor = corner_id>>latt_dir & 1; vector_incr(coords, this->lattice->getLatticeVector(latt_dir), factor * this->size * this->lattice->getConstant(latt_dir)); } Real distance_from_border = -this->domain->from_border(coords); if ( distance_from_border >= 0 ) { Real refinement = this->domain->interpolate(coords);// * this->lattice->getRepresentativeConstant(); /* got to satisfy two conditions to be split: * 1: not too close to border * 2: not larger than local refinement */ - bool cond1 = distance_from_border >= refinement; - bool cond2 = this->size > refinement; - if (cond1 && cond2) { + bool condition = this->size > refinement; + if (condition) { return true; } } } if (this->level <= 1) { return true; } return false; } template void Tree::create_points(PointContainer & points) { PointContainer temp_points; if (this->size == 0){ // in this case, the lattice needs to be inserted this->lattice->fillAtoms(this->south_west, temp_points); } else { //in this case, we have to set the south_west node and maybe neighbours temp_points.addPoint(south_west); if (this->neighbours != 0 && false) { Real neighbour_coords[DIM]; for (Uint neighbour = 0; neighbour < DIM ; ++neighbour) { if (!(this->neighbours>>neighbour & 1)) { vector_set(neighbour_coords, this->south_west); vector_incr(neighbour_coords, this->lattice->getLatticeVector(neighbour), this->size * this->lattice->getConstant(neighbour)); temp_points.addPoint(neighbour_coords); } } } } for (Uint i = 0; i < temp_points.getSize() ; ++i) { PointRef current_point = temp_points.getPoint(i); if (this->domain->inDomain(current_point)) { points.addPoint(current_point); } } } template void Tree::printself(std::ostream & stream, int indent) const { indent_space(stream, indent); stream << "IMPLEMENT THIS";} template unsigned int Tree::nb_boxes=0; template class Tree; template class Tree; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index a898a7c..32e05cb 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,44 +1,49 @@ #/** #* @file CMakeLists.txt #* @author Till Junge #* @date Thu Apr 12 15:05:28 2012 #* #* @brief #* #* @section LICENSE #* #* #* #*/ # #/* -------------------------------------------------------------------------- */ option(BLOCK_TEST "build the block test") if (BLOCK_TEST) register_test(block_test block_test.cc) endif(BLOCK_TEST) option(POINT_CONTAINER_TEST "build the point_container test") if (POINT_CONTAINER_TEST) register_test(point_container_test point_container_test.cc) endif(POINT_CONTAINER_TEST) option(FCC_LATTICE_TEST "build the fcc_lattice test") if (FCC_LATTICE_TEST) register_test(fcc_lattice_test fcc_lattice_test.cc) endif(FCC_LATTICE_TEST) option(MESH_TEST "build the mesh test") if (MESH_TEST) register_test(mesh_test mesh_test.cc) endif(MESH_TEST) option(MESH_TEST_2D "build the 2D mesh test") if (MESH_TEST_2D) register_test(mesh_test_2d mesh_test_2d.cc) endif(MESH_TEST_2D) option(DOMAIN_TEST "build the domain test") if (DOMAIN_TEST) register_test(domain_test domain_test.cc) endif(DOMAIN_TEST) + +option(DOMAIN_TEST_3D "build the 3D domain test") +if (DOMAIN_TEST_3D) + register_test(domain_test_3d domain_test_3d.cc) +endif(DOMAIN_TEST_3D) diff --git a/test/domain_test.cc b/test/domain_test.cc index f51c428..4dcf641 100644 --- a/test/domain_test.cc +++ b/test/domain_test.cc @@ -1,112 +1,112 @@ #include "cadd_mesher.hh" #include int main(int argc, char *argv[]) { const Uint dim = twoD; // Prepare a lattice //Real consts[dim] = {1, 1}; // Real consts[dim] = {4.04, 4.545}; Real consts[dim] = {5, 5}; //Real consts[dim] = {2.02, 2.23}; //Real miller[dim*dim] = {1, 1, -1, 1}; Real miller[dim*dim] = {1, 0, 0, 1}; FccLattice lattice(consts, miller); // Prepare the domain geometries Real overall_min_max[] = {-200, 200, -300, 0}; //Real overall_min_max[] = {-4.25, 4.25, -4.25, 4.25}; Block overall(overall_min_max, "overall"); Real refined_min_max[] = {-25, 25, -35, 0}; // Real refined_min_max[] = {-4.25, 4.25, -4.25, 4.25}; Block refined(refined_min_max, "refined"); Real atomistic_min_max[] = {-20, 20, -30, 0}; // Real atomistic_min_max[] = {-4.25, 4.25, -4.25, 4.25}; Block atomistic(atomistic_min_max, "atomistic"); - Real atomistic_skinned_min_max[] = {-16, 16, -26, 0}; + Real atomistic_skinned_min_max[] = {-19, 19, -29, 0}; // Real atomistic_skinned_min_max[] = {-16, 16, 26, 30}; Block atomistic_skinned(atomistic_skinned_min_max, "atomistic_skinned"); // Create domain Domain domain(&overall, &refined, &atomistic, &atomistic_skinned, &lattice); Real refinement_point[dim]; Real refinement; PointContainer punkte("FUER DEBUGGING"); std::vector aufloesung; // set up auxiliary mesh - for (Uint i = 0; i < dim ; ++i) { - for (Uint j = 0; j < dim ; ++j) { + for (int i = 0; i < dim ; ++i) { + for (int j = 0; j < dim ; ++j) { - refinement_point[0] = (1-2*static_cast(i))*200; - refinement_point[1] = -static_cast(j * 300); + refinement_point[0] = (1-2*static_cast(i))*201; + refinement_point[1] = -static_cast(j * 301 - (1-j)); refinement = 8; domain.setPointRefinement(refinement_point, refinement); punkte.addPoint(refinement_point); aufloesung.push_back(refinement); refinement_point[0] = (1-2*static_cast(i))*25; - refinement_point[1] = -static_cast(j * 35); + refinement_point[1] = -static_cast(j * 35 - (1-j)); refinement = 0; domain.setPointRefinement(refinement_point, refinement); refinement_point[0] = (1-2*static_cast(i))*100; - refinement_point[1] = -static_cast(j * 130); + refinement_point[1] = -static_cast(j * 130 - (1-j)); refinement = 8; domain.setPointRefinement(refinement_point, refinement); punkte.addPoint(refinement_point); aufloesung.push_back(refinement); } } std::cout << punkte << std::endl; Mesh testmesh(punkte.getVector()); testmesh.delaunay(1000); testmesh.initInterpolation(); std::cout << "Test mesh = " << std::endl << testmesh; // do all the domain test stuff here try { domain.fill_points(); } catch (std::string & error) { std::cout << error << std::endl; } std::cout << "number of tree boxes = " << domain.getTree().getNbBoxes() << std:: endl; domain.getPoints().dump("points_before_meshing"); domain.getAuxiliaryMesh().dump("domain_test_auxiliary_mesh"); Real point[][2] = {{-30, -5}, {-30, 0}, {-25, 0}, {-25, -5}}; for (Uint i = 0; i < 4 ; ++i) { Real * pt = point[i]; std::cout < int main(int argc, char *argv[]) { - const Uint dim = twoD; + const Uint dim = threeD; // Prepare a lattice //Real consts[dim] = {1, 1}; // Real consts[dim] = {4.04, 4.545}; - Real consts[dim] = {5, 5}; + Real consts[dim] = {5, 5, 5}; //Real consts[dim] = {2.02, 2.23}; //Real miller[dim*dim] = {1, 1, -1, 1}; - Real miller[dim*dim] = {1, 0, 0, 1}; + Real miller[dim*dim] = {1, 0, 0, 0, 1, 0, 0, 0, 1}; FccLattice lattice(consts, miller); // Prepare the domain geometries - Real overall_min_max[] = {-200, 200, -300, 0}; + Real overall_min_max[] = {-200, 200, -200, 200, -300, 0}; //Real overall_min_max[] = {-4.25, 4.25, -4.25, 4.25}; Block overall(overall_min_max, "overall"); - Real refined_min_max[] = {-25, 25, -35, 0}; + Real refined_min_max[] = {-25, 25, -25, 25, -35, 0}; // Real refined_min_max[] = {-4.25, 4.25, -4.25, 4.25}; Block refined(refined_min_max, "refined"); - Real atomistic_min_max[] = {-20, 20, -30, 0}; + Real atomistic_min_max[] = {-20, 20, -20, 20, -30, 0}; // Real atomistic_min_max[] = {-4.25, 4.25, -4.25, 4.25}; Block atomistic(atomistic_min_max, "atomistic"); - Real atomistic_skinned_min_max[] = {-16, 16, -26, 0}; + Real atomistic_skinned_min_max[] = {-19.9, 19.9, -19.9, 19.9, -29.9, 0}; // Real atomistic_skinned_min_max[] = {-16, 16, 26, 30}; Block atomistic_skinned(atomistic_skinned_min_max, "atomistic_skinned"); // Create domain Domain domain(&overall, &refined, &atomistic, &atomistic_skinned, &lattice); Real refinement_point[dim]; Real refinement; PointContainer punkte("FUER DEBUGGING"); std::vector aufloesung; // set up auxiliary mesh - for (Uint i = 0; i < dim ; ++i) { - for (Uint j = 0; j < dim ; ++j) { + for (int i = 0; i < dim-1; ++i) { + for (int j = 0; j < dim-1; ++j) { + for (int k = 0; k < dim-1; ++k) { - refinement_point[0] = (1-2*static_cast(i))*200; - refinement_point[1] = -static_cast(j * 300); - refinement = 8; - domain.setPointRefinement(refinement_point, refinement); - punkte.addPoint(refinement_point); - aufloesung.push_back(refinement); + refinement_point[0] = (1-2*static_cast(i))*201; + refinement_point[1] = (1-2*static_cast(j))*201; + refinement_point[2] = -static_cast(k * 301 - (1-k)); + refinement = 8; + domain.setPointRefinement(refinement_point, refinement); + punkte.addPoint(refinement_point); + aufloesung.push_back(refinement); - refinement_point[0] = (1-2*static_cast(i))*25; - refinement_point[1] = -static_cast(j * 35); - refinement = 0; + refinement_point[0] = (1-2*static_cast(i))*25; + refinement_point[1] = (1-2*static_cast(j))*25; + refinement_point[2] = -static_cast(k * 35 - (1-k)); + refinement = 0; - domain.setPointRefinement(refinement_point, refinement); + domain.setPointRefinement(refinement_point, refinement); - refinement_point[0] = (1-2*static_cast(i))*100; - refinement_point[1] = -static_cast(j * 130); - refinement = 8; + refinement_point[0] = (1-2*static_cast(i))*100; + refinement_point[1] = (1-2*static_cast(j))*100; + refinement_point[2] = -static_cast(k * 130 - (1-k)); + refinement = 8; - domain.setPointRefinement(refinement_point, refinement); - punkte.addPoint(refinement_point); - aufloesung.push_back(refinement); + domain.setPointRefinement(refinement_point, refinement); + punkte.addPoint(refinement_point); + aufloesung.push_back(refinement); + std::cout << "no of points = " << punkte.getSize() << std::endl; + } } } std::cout << punkte << std::endl; Mesh testmesh(punkte.getVector()); testmesh.delaunay(1000); testmesh.initInterpolation(); std::cout << "Test mesh = " << std::endl << testmesh; // do all the domain test stuff here try { domain.fill_points(); } catch (std::string & error) { std::cout << error << std::endl; } std::cout << "number of tree boxes = " << domain.getTree().getNbBoxes() << std:: endl; domain.getPoints().dump("points_before_meshing"); - domain.getAuxiliaryMesh().dump("domain_test_auxiliary_mesh"); - - Real point[][2] = {{-30, -5}, - {-30, 0}, - {-25, 0}, - {-25, -5}}; - for (Uint i = 0; i < 4 ; ++i) { - Real * pt = point[i]; - std::cout <