diff --git a/src/common/block_surface_manipulations.hh b/src/common/block_surface_manipulations.hh index 526daff..f3d7db0 100644 --- a/src/common/block_surface_manipulations.hh +++ b/src/common/block_surface_manipulations.hh @@ -1,219 +1,224 @@ #include template class Mesh; -template -inline Real distance2(const PointRef & pointA, const PointRef & pointB) { - Real dist2 = 0; - for (Uint i = 0; i < DIM ; ++i) { - Real manhattan = pointA.getComponent(i) - pointB.getComponent(i); - dist2 += manhattan * manhattan ; - } - return dist2; -} - -template -inline void add_if_ok(const PointRef & point, - const Mesh & auxiliary, - const std::vector & refinement, - const Real & step_size, - PointContainer & points) { - Real ref_sq = step_size * auxiliary.interpolate(point, - refinement); - Real step_sq = step_size*step_size; - ref_sq *= ref_sq; - if (ref_sq < step_sq) { - return; - } - for (Uint point_id = 0; point_id < points.getSize() ; ++point_id) { - if (ref_sq > distance2(point, points.getPoint(point_id))) { - return; - } - } - points.addPoint(point); -} - template void get_corners(const Real * min_max, PointContainer & corners, Uint direction = 0) { static Real point[DIM]; for (Uint i = 0; i < 2 ; ++i) { point[direction] = min_max[2*direction + i]; Uint new_dir = direction + 1; if (new_dir == DIM) { corners.addPoint(point); } else { get_corners(min_max, corners, new_dir); } } } template void get_edges(const PointContainer & corners, std::vector > & edge_ids) { Uint nb_corners = corners.getSize(); for (Uint i = 0; i < nb_corners ; ++i) { PointRef pointA = corners.getConstPoint(i); for (Uint j = i + 1 ; j < nb_corners ; ++j) { PointRefpointB = corners.getConstPoint(j); Uint diff = 0; for (Uint k = 0 ; k < DIM ; ++k) { diff += (pointA.getComponent(k) == pointB.getComponent(k)); } if (diff == DIM - 1) { edge_ids.push_back(std::pair(i, j)); } } } } template void get_surfaces(const PointContainer & corners, std::vector > & surface_ids) { throw(std::string("sorry, hardcoded")); } template<> void get_surfaces(const PointContainer & corners, std::vector > & surface_ids) { std::tuple surf; //surface x=xmin std::get<0>(surf) = 0; std::get<1>(surf) = 1; std::get<2>(surf) = 3; std::get<3>(surf) = 2; surface_ids.push_back(surf); //surface x=xmax std::get<0>(surf) = 4; std::get<1>(surf) = 5; std::get<2>(surf) = 7; std::get<3>(surf) = 6; surface_ids.push_back(surf); //surface y=ymin std::get<0>(surf) = 0; std::get<1>(surf) = 1; std::get<2>(surf) = 5; std::get<3>(surf) = 4; surface_ids.push_back(surf); //surface y=ymax std::get<0>(surf) = 2; std::get<1>(surf) = 3; std::get<2>(surf) = 7; std::get<3>(surf) = 6; surface_ids.push_back(surf); //surface z=zmin std::get<0>(surf) = 0; std::get<1>(surf) = 2; std::get<2>(surf) = 6; std::get<3>(surf) = 4; surface_ids.push_back(surf); //surface z=zmax std::get<0>(surf) = 1; std::get<1>(surf) = 3; std::get<2>(surf) = 7; std::get<3>(surf) = 5; surface_ids.push_back(surf); } template void add_corners(const PointContainer corners, const Mesh & auxiliary, const std::vector & refinement, const Real & step_size, PointContainer & points) { for (Uint i = 0; i < corners.getSize() ; ++i) { add_if_ok(corners.getConstPoint(i), auxiliary, refinement, step_size, points); } } template void add_edges(const PointContainer corners, const Mesh & auxiliary, const std::vector & refinement, const Real & step_size, PointContainer & points) { typedef std::vector > pair_vec; pair_vec edge_ids; get_edges(corners, edge_ids); pair_vec::iterator edge = edge_ids.begin(); pair_vec::iterator end_edge = edge_ids.end(); PointRef lower, upper; for (; edge != end_edge; ++edge) { PointRef point_a = corners.getConstPoint(edge->first); PointRef point_b = corners.getConstPoint(edge->second); Real length = pow(distance2(point_a, point_b), 0.5); Uint nb_steps = length/(2*step_size); Real step[DIM]; vector_difference(step, point_a.getArray(), point_b.getArray(), step_size/length); //these points will be iterated over the edge from both ends vector_set(lower.getArray(), point_a.getArray()); vector_set(upper.getArray(), point_b.getArray()); for (Uint i = 0 ; i < nb_steps ; ++i) { // move the lower bound by one step vector_incr(lower.getArray(), step); add_if_ok(lower, auxiliary, refinement, step_size, points); // move the upper bound by one step backwards vector_incr(upper.getArray(), step, -1.); add_if_ok(upper, auxiliary, refinement, step_size, points); } } } template void add_surfaces(const PointContainer & corners, const Mesh & auxiliary, const std::vector & refinement, const Real & step_size, PointContainer & points) { typedef std::tuple surf_type; std::vector surf_ids; get_surfaces(corners, surf_ids); std::vector::iterator surface = surf_ids.begin(); std::vector::iterator end_surface = surf_ids.end(); for (; surface != end_surface; ++surface) { PointContainer surf_corners; surf_corners.addPoint(corners.getConstPoint(std::get<0>(*surface))); surf_corners.addPoint(corners.getConstPoint(std::get<1>(*surface))); surf_corners.addPoint(corners.getConstPoint(std::get<2>(*surface))); surf_corners.addPoint(corners.getConstPoint(std::get<3>(*surface))); Real x_step[DIM], y_step[DIM]; Real x_length = pow(distance2(surf_corners.getPoint(0), surf_corners.getPoint(1)), 0.5); Real y_length = pow(distance2(surf_corners.getPoint(0), surf_corners.getPoint(3)), 0.5); vector_difference(x_step, surf_corners.getPoint(0).getArray(), surf_corners.getPoint(1).getArray(), step_size/x_length); vector_difference(y_step, surf_corners.getPoint(0).getArray(), surf_corners.getPoint(3).getArray(), step_size/y_length); Uint nb_steps = x_length/(2*step_size); keep_min(nb_steps, static_cast(y_length/(2*step_size))); for (Uint i = 0 ; i < nb_steps ; ++i) { vector_incr(surf_corners.getPoint(0).getArray(), x_step, 1); vector_incr(surf_corners.getPoint(0).getArray(), y_step, 1); vector_incr(surf_corners.getPoint(1).getArray(), x_step, -1); vector_incr(surf_corners.getPoint(1).getArray(), y_step, 1); vector_incr(surf_corners.getPoint(2).getArray(), x_step, -1); vector_incr(surf_corners.getPoint(2).getArray(), y_step, -1); vector_incr(surf_corners.getPoint(3).getArray(), x_step, 1); vector_incr(surf_corners.getPoint(3).getArray(), y_step, -1); add_corners(surf_corners, auxiliary, refinement, step_size, points); add_edges (surf_corners, auxiliary, refinement, step_size, points); } } } + + +template +void add_surface_points(const Real * min_max, + const Mesh & auxiliary, + const std::vector & refinement, + const Real step_size[DIM], + const bool periodicity[DIM], + PointContainer & points) { + PointRef point_hi, point_lo; + for (Uint xi = 0 ; xi < DIM ; ++xi) { + // We mesh from the outside inwards, to make sure that corners/edges get + // points + Uint ups = (xi+1) % DIM; + Real length_a = min_max[2*ups +1 ] - min_max[2*ups]; + Uint nb_steps_a = length_a/step_size[ups]; + for (Uint step_a = 0 ; step_a < nb_steps_a/2+1 ; ++step_a) { + point_lo[ups] = min_max[2*ups + 1] - step_a * step_size[ups]; + point_hi[ups] = min_max[2*ups] + step_a * step_size[ups]; + if (DIM == twoD) { + micro_add_surface_point(point_hi, point_lo, xi, min_max, auxiliary, + refinement, periodicity, points); + } else if (DIM == threeD) { + Uint zeta = (ups+1) % DIM; + Real length_b = min_max[2*zeta + 1] - min_max[2*zeta]; + Uint nb_steps_b = length_b/step_size[zeta]; + for (Uint step_b = 0 ; step_b < nb_steps_b/2+1 ; ++step_b) { + point_hi[zeta] = min_max[2*zeta + 1] - step_b * step_size[zeta]; + point_lo[zeta] = min_max[2*zeta] + step_b * step_size[zeta]; + micro_add_surface_point(point_hi, point_lo, xi, min_max, auxiliary, + refinement, periodicity, points); + } + } + } + } +} diff --git a/src/common/point_container.hh b/src/common/point_container.hh index 77f0e0f..7454e4e 100644 --- a/src/common/point_container.hh +++ b/src/common/point_container.hh @@ -1,379 +1,380 @@ /** * @file point_container.hh * @author Till Junge * @date Fri Apr 13 11:10:22 2012 * * @brief * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #ifndef __CADD_MESH_POINT_CONTAINER_HH__ #define __CADD_MESH_POINT_CONTAINER_HH__ #include #include "string" #include "sstream" #include "fstream" #include "common.hh" #include "dumper_paraview.hh" #include "dumper_lammps.hh" #include /*! \class PointRef simple representation of a point.*/ template class PointRef { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /*! Constructor * \param _vals array of which the point is part * \param _index position of the point in the array. E.g., for 3d points, the n-th point would have index n-1, and \a not 3*(n-1)*/ PointRef(const PointRef& other) :is_light(other.is_light) { if (this->is_light) { this->vals = other.vals; this->index = other.index; } else { this->vals = new Real [DIM]; std::copy(other.vals, other.vals + DIM, this->vals); this->index = 0; } }; PointRef & operator=(const PointRef& other) { if (this != &other) { if (this->is_light) { this->vals = new Real [DIM]; } std::copy(other.vals, other.vals + DIM, this->vals); this->is_light = false; this->index = 0; } return *this; } PointRef(Real * _vals=NULL, Uint _index=0) : index(_index) { if (_vals == NULL) { this->is_light = false; this->vals = new Real[DIM]; } else { this->is_light = true; this->vals = _vals + _index * DIM; } } PointRef(std::vector & coords) :index(0), is_light(false) { this->vals = new Real[DIM]; for (Uint i = 0; i < DIM ; ++i) { this->vals[i] = coords.at(i); } } bool operator==(const PointRef & other) { bool equal = true; for (Uint i = 0 ; i < DIM ; ++i) { equal &= (this->getComponent(i) == other.getComponent(i)); } return equal; } virtual ~PointRef(){ if (!this->is_light) delete[] this->vals;}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void make_me_independent() { if (this-> is_light) { Real new_vals[DIM]; for (Uint i = 0 ; i < DIM ; ++i) { new_vals[i] = this->vals[i]; } this->vals = new_vals; this->is_light = false; } } /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const { indent_space(stream, indent); stream << "Point "; stream.width(6); stream << this->index << ": "; stream << std::scientific << "("; for (Uint i = 0; i < DIM -1; ++i) { stream << this->vals[i] << ", "; } stream << this->vals[DIM -1] << ")" << std::endl; } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /*! self-explanatory */ + inline Real & operator[](Uint index) {return this->vals[index];} inline const Real & getComponent(Uint index) const {return this->vals[index];} inline Real * getArray() {return this->vals;} inline const Real * getArray() const {return this->vals;} inline std::vector getVector() const { return std::vector(this->vals, this->vals+DIM); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: Real * vals; Uint index; bool is_light; }; template class Geometry; /*! \class PointContainer container for vectors */ template class PointContainer { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /// Constructor \param name by default "" PointContainer(std::string _name = ""):name(_name){}; PointContainer(const PointContainer & other) :name(other.name), values(other.values) {}; virtual ~PointContainer(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: PointContainer & operator=(const PointContainer & other) { this->name = other.name; this->values = other.values; return *this; } /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const { indent_space(stream, indent); stream << "PointContainer '" << this->name << "' contains " << this->getSize() << " points:" << std::endl; Uint size = this->getSize(); Uint max_size = 20; Uint i; for (i = 0; i < size && i < max_size; ++i) { PointRef printref = this->getConstPoint(i); printref.printself(stream, indent + indent_incr); } if (size > max_size) { indent_space(stream, indent + indent_incr); stream << "Point ..." << std::endl; } } void computeBounds(Real * bounds) { for (Uint dir = 0 ; dir < DIM ; ++dir) { bounds[2 * dir + 0] = std::numeric_limits::max(); bounds[2 * dir + 1] = -std::numeric_limits::max(); } for (Uint point = 0 ; point < this->values.size() ; point+=DIM) { for (Uint dir = 0 ; dir < DIM ; ++dir) { keep_min(bounds[2*dir + 0], this->values.at(point+dir)); keep_max(bounds[2*dir + 1], this->values.at(point+dir)); } } } std::vector computeBounds() { std::vector bounds; bounds.resize(2*DIM); this->computeBounds(&bounds[0]); return bounds; } void filter(const Geometry & geom) { PointContainer temp; for (Uint point = 0 ; point < this->getSize() ; ++point) { PointRef current_point = this->getPoint(point); if (geom.is_inside(current_point)) { temp.addPoint(current_point); } } this->values = std::move(temp.values); } void dump_paraview(std::string filename="diagnostic_pointcontainer") const { iohelper::DumperParaview dumper; dumper.setMode(iohelper::TEXT); dumper.setPoints(const_cast(&this->values[0]), DIM, this->values.size()/DIM, filename); dumper.setConnectivity(NULL, iohelper::POINT_SET, 0, iohelper::C_MODE); dumper.setPrefix("./"); dumper.init(); try { dumper.dump(); } catch (iohelper::IOHelperException & err) { std::cout << err.what() << std::endl; throw(err); } } void dump_lammps(std::string filename, std::vector bounds) { this->dump_lammps(filename, &bounds[0]); } void dump_lammps(std::string filename="diagnostic_pointcontainer", Real * set_bounds = NULL) { Real bounds[2*DIM]; if (set_bounds == NULL) { this->computeBounds(bounds); } else { for (Uint i = 0 ; i < DIM*2 ; ++i) { bounds[i] = set_bounds[i]; } } iohelper::DumperLammps dumper(bounds); dumper.setPoints(const_cast(&this->values[0]), DIM, this->values.size()/DIM, filename); dumper.setPrefix("./"); dumper.setConnectivity(NULL, iohelper::POINT_SET, 0, iohelper::C_MODE); dumper.init(); dumper.dump(); } void dump_text(std::string filename="diagnostic_pointcontainer", std::string commentary="") const { std::ofstream ofile(filename.c_str()); if (ofile.good()) { ofile << "# file: " << filename <0) { ofile << "# " << commentary << std::endl; } } ofile << std::setprecision(16) << std::scientific << std::showpos; Uint current_point = 0; for (Uint i = 0 ; i < this->getSize(); ++i) { for (Uint j = 0 ; j < DIM - 1 ; ++j, ++current_point) { ofile << this->values.at(current_point) << "\t"; } ofile << this->values.at(current_point++); ofile << std::endl; } ofile.close(); } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /*! \param coords array of coordinates to add */ void addPoint(const Real * coords) { for (Uint i = 0; i < DIM; ++i) { this->values.push_back(coords[i]); } } void addPoint(std::vector & coords) throw (std::string) { if (coords.size() < DIM) { std::stringstream error_msg; error_msg << "Size of coords should be " << DIM << " but it's " << coords.size(); throw error_msg.str(); } this->addPoint(&coords[0]); } inline void addPoint(const PointRef & point) { this->addPoint(&point.getComponent(0));} /// returns the number of points in the container inline Uint getSize() const {return this->values.size()/DIM;} /// returns the container name std::string getName() const {return this->name;} /// sets the container name void setName(std::string _name) {this->name = _name;} /*! returns a PointRef to the point at position \a index. Throws a string if * the index is out of bounds * \param index */ PointRef getPoint(Uint index) throw (std::string) { if (index > this->getSize()) { std::stringstream error_msg; error_msg << "index " << index << " is out of bounds. Size of vector " << this->name << " is only " << this->getSize(); throw error_msg.str(); } return PointRef(&this->values[0], index); } const PointRef getConstPoint(Uint index) const throw (std::string) { if (index > this->getSize()) { std::stringstream error_msg; error_msg << "index " << index << " is out of bounds. Size of vector " << this->name << " is only " << this->getSize(); throw error_msg.str(); } const PointRef ret_ref(const_cast(&values[0]), index); return ret_ref; } /*! returnd the vector of raw data. You'll probably never use this.*/ const std::vector & getVector() const {return this-> values;} /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: std::string name; std::vector values; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator <<(std::ostream & stream, const PointContainer & _this) { _this.printself(stream); return stream; } /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator <<(std::ostream & stream, const PointRef & _this) { _this.printself(stream); return stream; } #endif /* __CADD_MESH_POINT_CONTAINER_HH__ */ diff --git a/src/domain/domain.cc b/src/domain/domain.cc index 82ee698..25eb8af 100644 --- a/src/domain/domain.cc +++ b/src/domain/domain.cc @@ -1,251 +1,252 @@ #include "domain.hh" #include "tree.hh" #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) :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()) { 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() throw(std::string){ +void Domain::fill_points(bool periodicity[DIM]) throw(std::string){ this->initialiseAuxiliaryMesh(); if (this->boundaries_special_treatment) { this->overall->generate_surface_points(*this->auxiliary_mesh, this->refinement, - this->lattice->getRepresentativeConstant(), + this->lattice->getConstants(), + periodicity, this->points); } this->tree->fill(this->points); this->fill_atoms(); this->filled = true; } 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() { this->auxiliary_mesh = new Mesh(this->refinement_points.getVector()); this->auxiliary_mesh->delaunay(1000); 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->atomistic_skinned); try { this->main_mesh->delaunay(quality); } catch (std::string & error) { this->tree->dump(); this->points.dump_paraview(); throw("affen_fehler " + error); } this->main_mesh->cleanup(*this->atomistic); this->main_meshed = true; } return *this->main_mesh; } template 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() 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 = int(this->interpolate(point)) * this->lattice->getRepresentativeConstant(); bool retval; if (this->boundaries_special_treatment) { retval = distance_from_border >= refinement && !this->atomistic_skinned->is_inside(point); } else { retval = !this->atomistic_skinned->is_inside(point); } return retval; } template void Domain::fill_atoms() { Uint current_index = 0; 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->atoms.addPoint(point); if (this->inFemDomain(point)) { this->pad_atoms_indices.push_front(current_index); } ++ current_index; } } } template void Domain::dump_atoms_lammps(std::string filename) { if (this->atoms.getSize() == 0) { this->fill_atoms(); } this->atoms.dump_lammps(filename); } template void Domain::dump_atoms_paraview(std::string filename) { if (this->atoms.getSize() == 0) { this->fill_atoms(); } this->atoms.dump_paraview(filename); } template const PointContainer & Domain::getAtoms() { if (this->atoms.getSize() == 0) { this->fill_atoms(); } return this->atoms; } 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) { int int_refinement = this->interpolate(point); Real refinement = int_refinement * this->lattice->getRepresentativeConstant(); return distance_from_border >= 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) { this->fill_atoms(); } if (&atomistic == NULL) { this->getMesh().compute_interface_atoms(interface, *this->refined); } else { this->getMesh().compute_interface_atoms(interface, atomistic); } // now, clean the pad atoms to get rid of the interface atoms for (Uint interface_index = 0 ; interface_index < interface.getSize() ; ++interface_index) { PointRef interface_point = interface.getPoint(interface_index); auto pad_index = this->pad_atoms_indices.begin(); auto pad_before_index = this->pad_atoms_indices.before_begin(); auto pad_end = this->pad_atoms_indices.end(); bool searching = true; for (; pad_index != pad_end && searching; ++pad_index, ++pad_before_index) { PointRef pad_point = this->atoms.getPoint(*pad_index); if (pad_point == interface_point) { this->pad_atoms_indices.erase_after(pad_before_index); searching=false; } } } auto pad_index = this->pad_atoms_indices.begin(); auto pad_end = this->pad_atoms_indices.end(); for(;pad_index != pad_end; ++pad_index) { pad.addPoint(this->atoms.getPoint(*pad_index)); } } template class Domain; template class Domain; diff --git a/src/domain/domain.hh b/src/domain/domain.hh index 0e90265..6ec2fc8 100644 --- a/src/domain/domain.hh +++ b/src/domain/domain.hh @@ -1,212 +1,222 @@ /** * @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- * 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(const Geometry & overall_, const Geometry & refined_, const Geometry & atomistic_, const Geometry & atomistic_skinned_, const Lattice & lattice_, bool boundaries_special_treatment = 0, const PointRef & centre=*reinterpret_cast *>(NULL)); 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() throw(std::string); + 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(const std::vector & periodicity) throw(std::string) { + bool per[DIM] = {0}; + for (Uint i = 0 ; i < DIM ; ++i) { + per[i] = periodicity.at(0); + } + return this->fill_points(per); + } 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(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)));} /*! 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);} /* ------------------------------------------------------------------------ */ /* 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() throw(std::string); 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(); // 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 filles when the atoms are /// computed, so that we can separate interface and pad std::forward_list pad_atoms_indices; /// 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; }; /* -------------------------------------------------------------------------- */ /* 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/geometry/block.cc b/src/geometry/block.cc index 13b306a..87d533b 100644 --- a/src/geometry/block.cc +++ b/src/geometry/block.cc @@ -1,230 +1,346 @@ #include "block.hh" #include "block_surface_manipulations.hh" #include #include #include #include #include #include #include #include #include class BlockException: public std::exception { public: BlockException(std::stringstream & _what):what_stream(_what.str()) {}; BlockException(std::string &_what):what_stream(_what) {}; ~BlockException() throw() {}; virtual const char* what() const throw() { return this->what_stream.c_str(); } private: std::string what_stream; }; template Block::Block(const std::vector & _min_max, std::string _name) :Geometry(_name) { if (_min_max.size() != 2*DIM) { std::stringstream error; error << "The length of 'min_max' should be = " << 2*DIM; throw BlockException(error); } for (Uint i = 0 ; i < 2*DIM ; ++i) { this->min_max[i] = _min_max[i]; } } template Block::Block(const Real * min_max, std::string _name) :Geometry(_name) { for (Uint i = 0; i < 2 * DIM; ++i) { this->min_max[i] = min_max[i]; } } template Block::Block(const Block & other) :Geometry(other.name) { for (Uint i = 0 ; i < 2 * DIM ; ++i) { this->min_max[i] = other.min_max[i]; } } template void Block::shift(const PointRef & offset) { for (Uint i = 0 ; i < DIM ; ++i) { this->min_max[2 * i + 0] += offset.getComponent(i); this->min_max[2 * i + 1] += offset.getComponent(i); } } template Real Block::from_border(const Real * point) const { Real max_val = -std::numeric_limits::max(); for (Uint i = 0; i < DIM; ++i) { keep_max(max_val, this->min_max[2*i] - point[i]); keep_max(max_val, point[i] - this->min_max[2*i+1]); } return max_val; } template inline void in_direction_helper(Real & min_val, Real & max_val, const Real * direction, Real * point) { Real dot_prod = 0; Real norm = 0; for (Uint i = 0; i < DIM; ++i) { dot_prod += point[i] * direction[i]; norm += direction[i] * direction[i]; } dot_prod /= sqrt(norm); if (eval_min) { keep_min(min_val, dot_prod); } if (eval_max) { keep_max(max_val, dot_prod); } } template Real Block::min_in_direction(const Real * dir, const Real & unit) const { Real point[DIM]; Real min_val = std::numeric_limits::max(); for (Uint i = 0; i < 2; ++i) { point[0] = min_max[i]; for (Uint j = 0; j < 2; ++j) { point[1] = min_max[2 + j]; if (DIM == 3) { for (Uint k = 0; k < 2; ++k) { point[2] = min_max[4 + k]; in_direction_helper(min_val, min_val, dir, point); } } else if (DIM == 2) { in_direction_helper(min_val, min_val, dir, point); } } } return min_val/unit; } template Real Block::max_in_direction(const Real * dir, const Real & unit) const { Real point[DIM]; Real max_val = -std::numeric_limits::max(); for (Uint i = 0; i < 2; ++i) { point[0] = min_max[i]; for (Uint j = 0; j < 2; ++j) { point[1] = min_max[2 + j]; if (DIM == 3) { for (Uint k = 0; k < 2; ++k) { point[2] = min_max[4 + k]; in_direction_helper(max_val, max_val, dir, point); } } else if (DIM == 2) { in_direction_helper(max_val, max_val, dir, point); } } } return max_val/unit; } template Real Block::size_in_direction(const Real * dir, const Real & unit) const { Real point[DIM]; Real min_val = std::numeric_limits::max(); Real max_val = -std::numeric_limits::max(); for (Uint i = 0; i < 2; ++i) { point[0] = min_max[i]; for (Uint j = 0; j < 2; ++j) { point[1] = min_max[2 + j]; if (DIM == 3) { for (Uint k = 0; k < 2; ++k) { point[2] = min_max[4 + k]; in_direction_helper(min_val, max_val, dir, point); } } else if (DIM == 2) { in_direction_helper(min_val, max_val, dir, point); } } } return (max_val - min_val)/unit; } template void Block::printself(std::ostream & stream, int indent) const { indent_space(stream, indent); stream << "Block " << this->get_name() << ": " << std::endl; indent_space(stream, indent); stream << "x_min " << this->min_max[0] << ", x_max = " << this->min_max[1] << std::endl; indent_space(stream, indent); stream << "y_min " << this->min_max[2] << ", y_max = " << this->min_max[3] << std::endl; if (DIM == 3) { indent_space(stream, indent); stream << "z_min " << this->min_max[4] << ", z_max = " << this->min_max[5] << std::endl; } } -template<> -void Block::generate_surface_points(const Mesh & auxiliary, - const std::vector & refinement, - const Real & step_size, - PointContainer & points) const { - PointContainer corners; - get_corners(this->min_max, corners); - add_corners(corners, auxiliary, refinement, step_size, points); - add_edges (corners, auxiliary, refinement, step_size, points); -} - -template<> -void Block::generate_surface_points(const Mesh & auxiliary, - const std::vector & refinement, - const Real & step_size, - PointContainer & points) const { - PointContainer corners; - get_corners(this->min_max, corners); - add_corners (corners, auxiliary, refinement, step_size, points); - add_edges (corners, auxiliary, refinement, step_size, points); - add_surfaces(corners, auxiliary, refinement, step_size, points); +//template<> +//void Block::generate_surface_points(const Mesh & auxiliary, +// const std::vector & refinement, +// const Real & step_size, +// PointContainer & points, +// bool periodicity[twoD]) const { +// PointContainer corners; +// get_corners(this->min_max, corners); +// add_corners(corners, auxiliary, refinement, step_size, points); +// add_edges (corners, auxiliary, refinement, step_size, points, periodicity); +//} +// +//template<> +//void Block::generate_surface_points(const Mesh & auxiliary, +// const std::vector & refinement, +// const Real & step_size, +// PointContainer & points, +// bool periodicity[threeD]) const { +// PointContainer corners; +// get_corners(this->min_max, corners); +// add_corners (corners, auxiliary, refinement, step_size, points); +// add_edges (corners, auxiliary, refinement, step_size, points, periodicity); +// add_surfaces(corners, auxiliary, refinement, step_size, points, periodicity); +//} + +template +inline Real distance2(const PointRef & pointA, const PointRef & pointB) { + Real dist2 = 0; + for (Uint i = 0; i < DIM ; ++i) { + Real manhattan = pointA.getComponent(i) - pointB.getComponent(i); + dist2 += manhattan * manhattan ; + } + return dist2; +} + +template +inline bool add_if_ok(const PointRef & point, + const Mesh & auxiliary, + const std::vector & refinement, + const Real & step_size, + PointContainer & points) { + Real ref_sq = step_size * auxiliary.interpolate(point, + refinement); + Real step_sq = step_size*step_size; + ref_sq *= ref_sq; + if (ref_sq < step_sq) { + return false; + } + for (Uint point_id = 0; point_id < points.getSize() ; ++point_id) { + if (ref_sq > distance2(point, points.getPoint(point_id))) { + return false; + } + } + points.addPoint(point); + return true; +} + +template +inline void micro_add_surface_point(PointRef & point_hi, + PointRef & point_lo, + const Uint xi, + const Real * min_max, + const Mesh & auxiliary, + const std::vector & refinement, + const bool periodicity[DIM], + const Real step_size[DIM], + PointContainer & points) { + point_hi[xi] = point_lo[xi] = min_max[2*xi]; + bool hi_ok = add_if_ok(point_hi, auxiliary, refinement, step_size[DIM], points); + bool lo_ok = add_if_ok(point_lo, auxiliary, refinement, step_size[DIM], points); + point_hi[xi] = point_lo[xi] = min_max[2*xi+1]; + if (!periodicity[xi]) { + add_if_ok(point_hi, auxiliary, refinement, step_size[DIM], points); + add_if_ok(point_lo, auxiliary, refinement, step_size[DIM], points); + } else { + if (hi_ok) points.addPoint(point_hi); + if (lo_ok) points.addPoint(point_lo); + } } +template +void recurse_through_dimensions(PointRef & point_hi, + PointRef & point_lo, + const Uint xi, + const Uint my_offset_dim, + const Real * min_max, + const Mesh & auxiliary, + const std::vector & refinement, + const bool periodicity[DIM], + const Real step_size[DIM], + PointContainer & points, + Real * starts) { + Uint my_dim = (xi+my_offset_dim) % DIM; + Real length = min_max[2*my_dim + 1] - min_max[2*my_dim]; + Uint nb_steps = length/step_size[my_dim]; + starts[DIM-1-xi] = 1; + for (Uint step = starts[my_offset_dim-1]; step < nb_steps/2+1; ++step) { + + } +} + +template +void Block::generate_surface_points(const Mesh & auxiliary, + const std::vector & refinement, + const Real step_size[DIM], + const bool periodicity[DIM], + PointContainer & points) const { + PointRef point_hi, point_lo; + Uint starts[DIM] = {0}; + + for (Uint xi = 0 ; xi < DIM ; ++xi) { + // We mesh from the outside inwards, to make sure that corners/edges get + // points + Uint ups = (xi+1) % DIM; + Real length_a = this->min_max[2*ups +1 ] - this->min_max[2*ups]; + Uint nb_steps_a = length_a/step_size[ups]; + starts[DIM-1-xi] = 1; + for (Uint step_a = starts[0] ; step_a < nb_steps_a/2+1 ; ++step_a) { + point_lo[ups] = this->min_max[2*ups + 1] - step_a * step_size[ups]; + point_hi[ups] = this->min_max[2*ups] + step_a * step_size[ups]; + if (DIM == twoD) { + micro_add_surface_point(point_hi, point_lo, xi, this->min_max, auxiliary, + refinement, periodicity, step_size, points); + } else if (DIM == threeD) { + Uint zeta = (ups+1) % DIM; + Real length_b = this->min_max[2*zeta + 1] - this->min_max[2*zeta]; + Uint nb_steps_b = length_b/step_size[zeta]; + for (Uint step_b = starts[1] ; step_b < nb_steps_b/2+1 ; ++step_b) { + point_hi[zeta] = this->min_max[2*zeta + 1] - step_b * step_size[zeta]; + point_lo[zeta] = this->min_max[2*zeta] + step_b * step_size[zeta]; + micro_add_surface_point(point_hi, point_lo, xi, this->min_max, auxiliary, + refinement, periodicity, step_size, points); + } + } + } + } +} + + template void Block::complement_periodically(const Mesh & auxiliary, const std::vector & refinement, PointContainer & points, int direction, Uint dim) const { // relevant coordinate of the good and the bad side Real ok_coord, bad_coord; if (direction > 0) { ok_coord = this->min_max[2*dim + 0]; bad_coord = this->min_max[2*dim + 1]; } else { ok_coord = this->min_max[2*dim + 1]; bad_coord = this->min_max[2*dim + 0]; } Real tol = std::max(fabs(ok_coord), 1.); tol = 2*tol*std::numeric_limits::epsilon(); PointRef point; for (Uint i = 0; i < points.getSize() ; ++i) { point = points.getPoint(i); if (fabs(point.getComponent(dim) - ok_coord) < tol) { point.getArray()[dim] = bad_coord; Real local_refinement = auxiliary.interpolate(point, refinement); if (local_refinement > 0) { points.addPoint(point); } } } } template class Block<2>; template class Block<3>; diff --git a/src/geometry/block.hh b/src/geometry/block.hh index 490d4d8..1d1a9b6 100644 --- a/src/geometry/block.hh +++ b/src/geometry/block.hh @@ -1,105 +1,106 @@ /** * @file block.hh * @author Till Junge * @date Wed Apr 4 15:28:40 2012 * * @brief Block geometry class * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #ifndef __cadd_mesh_BLOCK_HH__ #define __cadd_mesh_BLOCK_HH__ #include "common.hh" #include "geometry.hh" /*! \brief implementation of the Geometry interface for blocks aligned with coordinate axes*/ template class Block:public Geometry { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: // min_max is x_min, x_max, y_min, y_max ... Block(const std::vector & min_max, std::string name = ""); Block(const Real * min_max, std::string _name = ""); Block(const Block & other); virtual ~Block(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual Geometry * resolveType() const { return newGeom(*this); } /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; virtual void shift (const PointRef & offset); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ /* Returns the minimum value of all constraint violations (negative return value means point is inside and return value from the nearest boundary) */ virtual bool is_inside(const Real * point) const { return this->from_border(point) <= 0;} virtual Real from_border(const Real * point) const; // for instance x_min of a geom in latticespacings virtual Real min_in_direction(const Real * dir, const Real & unit = 1) const; // idem virtual Real max_in_direction(const Real * dir, const Real & unit = 1) const; // idem virtual Real size_in_direction(const Real * dir, const Real & unit = 1) const; virtual void generate_surface_points(const Mesh & auxiliary, const std::vector & refinement, - const Real & step_size, + const Real step_size[DIM], + const bool periodicity[DIM], PointContainer & points) const; virtual void complement_periodically(const Mesh & auxiliary, const std::vector & refinement, PointContainer & points, int direction, Uint dim) const; std::vector bounds () { std::vector bounds; for (Uint i = 0 ; i < DIM*2 ; ++i) { bounds.push_back(this->min_max[i]); } return bounds; } private: Real min_max[2*DIM]; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator <<(std::ostream & stream, const Block & _this) { _this.printself(stream); return stream; } #endif /* __cadd_mesh_BLOCK_HH__ */ diff --git a/src/geometry/geometry.hh b/src/geometry/geometry.hh index e3b75d9..0f0be3e 100644 --- a/src/geometry/geometry.hh +++ b/src/geometry/geometry.hh @@ -1,132 +1,133 @@ /** * @file geometry.hh * @author Till Junge * @date Wed Apr 4 15:15:10 2012 * * @brief Geometry motherclass (interface) for CADD mesher * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #ifndef __GEOMETRY_HH__ #define __GEOMETRY_HH__ #include "common.hh" #include "point_container.hh" #include //#include "cadd_mesh.hh" template class Mesh; /*! \brief Geometry interface for CADD mesher * You will never use this class.*/ template class Geometry { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /*! Constructor *\param _name a string identifier for the geometry */ Geometry(std::string _name = ""):name(_name){}; virtual ~Geometry(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// function to print the content of the class virtual void printself(std::ostream & stream, int indent = 0) const { stream << "This is just a interface class and this function should never " << "ever be called" < & offset) = 0; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ /// returns the string identifier const std::string & get_name() const { return this->name; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ public: virtual Geometry * resolveType() const = 0; /*! test whether a point is inside the domain * \param point coordinates of the point to test*/ virtual bool is_inside(const Real * point) const = 0; inline bool is_inside(const PointRef & point) const { return this->is_inside(&(point.getComponent(0)));} /*! Returns the minimum value of all constraint violations (negative return * value means point is inside and return value from the nearest boundary) * \param point coordintates of the point to test*/ virtual Real from_border(const Real * point) const = 0; inline Real from_border(const PointRef & point) const { return this->from_border(&(point.getComponent(0)));} /*! projects the geometry on a line of direction \a dir and returns the * minimum of the line segment spanned * \param dir direction vector (will be scaled to unity by the method * \param unit unit in which the distance should be expressed, by default unity*/ virtual Real min_in_direction(const Real * dir, const Real & unit = 1) const = 0; /*! projects the geometry on a line of direction \a dir and returns the * maximum of the line segment spanned * \param dir direction vector (will be scaled to unity by the method * \param unit unit in which the distance should be expressed, by default unity*/ virtual Real max_in_direction(const Real * dir, const Real & unit = 1) const = 0; /*! projects the geometry on a line of direction \a dir and returns the length * of the projection * \param dir direction vector (will be scaled to unity by the method * \param unit unit in which the distance should be expressed, by default unity*/ virtual Real size_in_direction(const Real * dir, const Real & unit = 1) const = 0; virtual void generate_surface_points(const Mesh & auxiliary, const std::vector & refinement, - const Real & step_size, + const Real step_size[DIM], + const bool periodicity[DIM], PointContainer & points) const = 0; /*! When the dimensions of the overall geometry are not a power of two times * the lattice size, the domain can be not filled correctly. This method * corrects the fucked up side by periodically repeating the ok side. * \param direction positive values mean that the min side is used to augment * the max size, negative values the inverse. * \param dim dimension to fix: 0:x, 1:y, 2:z */ virtual void complement_periodically(const Mesh & auxiliary, const std::vector & refinement, PointContainer & points, int direction, Uint dim) const = 0; protected: std::string name; }; /// standard output stream operator template inline std::ostream & operator <<(std::ostream & stream, const Geometry & _this) { _this.printself(stream); return stream; } template Geometry * newGeom(const subGeom & geometry) { return static_cast*>(new subGeom(geometry)); } #endif /* __GEOMETRY_HH__ */ diff --git a/src/swig/CaddMesher.i b/src/swig/CaddMesher.i index 3730e63..857daea 100644 --- a/src/swig/CaddMesher.i +++ b/src/swig/CaddMesher.i @@ -1,160 +1,160 @@ %module CADDMesher - +#pragma SWIG nowarn=325,389 %{ #include "../src/common/cadd_mesher.hh" %} typedef double Real; typedef unsigned int Uint; %include "std_string.i" %include "std_vector.i" %include "typemaps.i" %define __STR__() std::string __str__() { std::ostringstream out; out << *$self; return out.str(); } %enddef namespace std { %template(vecReal) std::vector ; %template(vecInt) std::vector ; + %template(vecBool) std::vector ; } %rename(__assign__) *::operator=; %apply Real &INPUT {Real & refinement}; %include "../src/domain/domain.hh" //-------------------------------------------- // Domain %extend Domain { __STR__()}; %ignore Domain::setPointRefinement(Real const *, Real &); %template(Domain2d) Domain; %extend Domain { __STR__()}; %ignore Domain::setPointRefinement(Real const *, Real &); %template(Domain3d) Domain; //-------------------------------------------- // Geometry %include "../src/geometry/geometry.hh" %extend Geometry { __STR__() }; %template(Geometry2d) Geometry; %extend Geometry { __STR__() }; %template(Geometry3d) Geometry; %include "../src/geometry/block.hh" %ignore Block::Block( Real const *, std::string); %template(Block2d) Block; %ignore Block::Block( Real const *, std::string); %template(Block3d) Block; //-------------------------------------------- // Lattice %include "../src/lattice/lattice.hh" %extend Lattice { __STR__() }; %template(Lattice2d) Lattice; %extend Lattice { __STR__() }; %template(Lattice3d) Lattice; %include "../src/lattice/fcc_lattice.hh" %ignore FccLattice::FccLattice(Real *, Real *); %ignore FccLattice::FccLattice(Real *); %template(FccLattice2d) FccLattice; %ignore FccLattice::FccLattice(Real *, Real *); %ignore FccLattice::FccLattice(Real *); %template(FccLattice3d) FccLattice; //-------------------------------------------- // PointContainer and PointRefs %include "../src/common/point_container.hh" %ignore PointContainer::addPoint(const Real*); %extend PointContainer { __STR__() }; %template(PointContainer2d) PointContainer; %ignore PointContainer::addPoint(const Real*); %extend PointContainer { __STR__() }; %template(PointContainer3d) PointContainer; %ignore PointRef::operator=; %ignore PointRef::PointRef(Real *, Uint); %ignore PointRef::PointRef(Real *); %extend PointRef { __STR__()}; %extend PointRef { inline void __setitem__(size_t i, const Real & val) throw(std::out_of_range) { if (i>=twoD || i<0) { throw std::out_of_range("out of bounds access"); } $self->getArray()[i] = val; } } %extend PointRef { inline double __getitem__(size_t i) throw(std::out_of_range) { if (i>=twoD || i<0) { throw std::out_of_range("out of bounds access"); } return $self->getArray()[i]; } } %template (PointRef2d) PointRef; - - %ignore PointRef::PointRef(Real *, Uint); %ignore PointRef::PointRef(Real *); %extend PointRef { __STR__()}; %extend PointRef { inline void __setitem__(size_t i, const Real & val) throw(std::out_of_range) { if (i>=threeD || i<0) { throw std::out_of_range("out of bounds access"); } $self->getArray()[i] = val; } } %extend PointRef { inline double __getitem__(size_t i) throw(std::out_of_range) { if (i>=threeD || i<0) { throw std::out_of_range("out of bounds access"); } return $self->getArray()[i]; } } %template (PointRef3d) PointRef; //-------------------------------------------- // Mesh %include "../src/domain/cadd_mesh.hh" %extend Mesh { __STR__()}; %template (Mesh2d) Mesh; %extend Mesh { __STR__()}; %template (Mesh3d) Mesh; //-------------------------------------------- // Tree + %include "../src/domain/tree.hh" %extend Tree { __STR__()}; %ignore Tree::iterator; %ignore Tree::getIterator(); %template (Tree2d) Tree; %extend Tree { __STR__()}; %ignore Tree::iterator; %ignore Tree::getIterator(); %template (Tree3d) Tree; %include "../src/common/cadd_mesher.hh"