diff --git a/src/dd/ref_elem_dd.hh b/src/dd/ref_elem_dd.hh index 57101d2..63b4f41 100644 --- a/src/dd/ref_elem_dd.hh +++ b/src/dd/ref_elem_dd.hh @@ -1,138 +1,148 @@ /** * @file ref_node_dd.hh * * @author Guillaume Anciaux * @author Till Junge * @author Max Hodapp * @author Jaehyun Cho * * @date Fri Nov 15 14:49:04 2013 * * @brief Generic reference to DD nodes * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * LibMultiScale is free software: you can redistribute it and/or modify it * under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * LibMultiScale is distributed in the hope that it will be useful, but * WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with LibMultiScale. If not, see . * */ #ifndef __LIBMULTISCALE_REF_ELEM_DD_HH__ #define __LIBMULTISCALE_REF_ELEM_DD_HH__ /* ----------------------------------------------------------------- */ #include "container_mesh.hh" #include "lm_common.hh" #include "ref_element.hh" #include "ref_node_dd.hh" /* ----------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* ----------------------------------------------------------------- */ template struct AccessorDDElemDof : public AccessorDof { using AccessorDof::AccessorDof; using AccessorDof::operator=; }; /* -------------------------------------------------------------------------- */ // template struct getField { // template static inline auto get(Ref &&ref) { // LM_TOIMPLEMENT; // return Vector<1>(); // } // }; // template <> struct getField<_burgers> { // template static inline decltype(auto) get(Ref &&ref) { // return ref.burgers(); // } // }; // template <> struct getField<_normal> { // template static inline decltype(auto) get(Ref &&ref) { // return ref.normal(); // } // }; /* ----------------------------------------------------------------- */ /** * Class RefDDElem * */ template class RefDDElem : public RefElement { public: static constexpr UInt Dim = _Dim; -#define declare_accessor(_name) \ +#define declare_accessor(_name) \ auto _name() { return static_cast(this)->_name(); } declare_accessor(burgers); declare_accessor(normal); #undef declare_accessor std::vector globalIndexes() override { LM_TOIMPLEMENT; }; template AccessorDDElemDof, ftype> field() { AccessorDDElemDof, ftype> acc(*this); return acc; } }; /* ------------------------------------------------------------------------ */ template class RefGenericDDElem : public RefDDElem<_Dim, RefGenericDDElem<_Dim>> { public: using fields = field_list<_burgers, _normal>; - - Vector<_Dim> burgers() { LM_TOIMPLEMENT; } - Vector<_Dim> normal() { LM_TOIMPLEMENT; } + + Vector<_Dim> &burgers() { return burg; } + Vector<_Dim> &normal() { return norm; } std::vector globalIndexes() override { LM_TOIMPLEMENT; } - std::vector localIndexes() override { LM_TOIMPLEMENT; } + std::vector localIndexes() override { return conn; } template void setConnectivity(Vec &&conn) { - auto n_nodes = conn.rows() * conn.cols(); this->conn.clear(); - for (UInt i = 0; i < n_nodes; ++i) - this->conn.push_back(conn[i]); + for (auto &&idx : conn) + this->conn.push_back(idx); + } + + template void setNodes(Vec &nodes) { + this->nodes.clear(); + for (auto &&idx : this->conn) { + this->nodes.push_back(&nodes[idx]); + } } - RefGenericDDNode<_Dim> &getNode(UInt) { LM_TOIMPLEMENT; } + RefGenericDDNode<_Dim> &getNode(UInt idx) { return *nodes[idx]; } bool operator==(RefGenericDDElem<_Dim> &) { LM_TOIMPLEMENT; } RefGenericDDElem getMirrorElem() { LM_TOIMPLEMENT; } +private: std::vector conn; + std::vector *> nodes; + Vector<_Dim> burg; + Vector<_Dim> norm; }; /* ------------------------------------------------------------------------ */ template using ContainerGenericDDMesh = ContainerMesh>, ContainerArray>>; __END_LIBMULTISCALE__ #endif /* __LIBMULTISCALE_REF_ELEM_DD_HH__ */ diff --git a/src/filter/compute_cadd_template.cc b/src/filter/compute_cadd_template.cc index a57ce0e..5e0b7a5 100644 --- a/src/filter/compute_cadd_template.cc +++ b/src/filter/compute_cadd_template.cc @@ -1,418 +1,441 @@ #include "compute_cadd_template.hh" #include "factory_multiscale.hh" #include "lib_dd.hh" #include "lib_dumper.hh" #include "lib_filter.hh" #include "lib_md.hh" #include "lib_stimulation.hh" #include "lm_common.hh" #include #include #include #include #include #include __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ void ComputeCADDTemplate::loadNPZ(std::string &filename) { py::object load = py::module::import("numpy").attr("load"); auto npz = load(filename); auto field_names = npz.attr("files").cast>(); auto fields = npz.cast().cast>(); // std::for_each(field_names.begin(), field_names.end(), // [](auto &&v) { std::cout << v << std::endl; }); auto triangles = fields["triangles"].cast>(); auto points = fields["points"].cast>(); auto corrections = fields["correction"].cast>(); auto &nodes = mesh.getContainerNodes(); RefPointData<2> p; for (UInt i = 0; i < points.rows(); ++i) { p = VectorView<2>(&points(i, 0)); nodes.push_back(p); } for (UInt i = 0; i < corrections.rows(); ++i) { auto corr = VectorView<3>(&corrections(i, 0)); this->corrections.push_back(corr); } auto &elems = mesh.getContainerElems(); for (UInt i = 0; i < triangles.rows(); ++i) { RefGenericElem<2> el; for (UInt n = 0; n < 3; ++n) el.conn.push_back(triangles(i, n)); elems.push_back(el); } } /* -------------------------------------------------------------------------- */ ComputeCADDTemplate::ComputeCADDTemplate(const LMID &name) : LMObject(name), mesh(this->getID() + ":mesh") { this->removeInput("input1"); this->createInput("points"); this->createInput("segments"); + this->createOutput("displacements") = + ContainerArray(this->getID() + ":" + "displacements"); + this->inversion = TOP_RIGHT; this->search_dim = 1; this->tol = -1; this->r_cut = -1; this->factor[0] = this->factor[1] = 1; } /* -------------------------------------------------------------------------- */ ComputeCADDTemplate::~ComputeCADDTemplate() {} /* -------------------------------------------------------------------------- */ inline bool is_in_triangle(Vector<2> &a, Vector<2> &b, Vector<2> &c, Vector<2> p) { Vector<2> v0; Vector<2> v1; Vector<2> v2; Real u, v; v0 = c - a; v1 = b - a; v2 = p - a; // Compute dot products Real dot00 = v0.dot(v0); Real dot01 = v0.dot(v1); Real dot02 = v0.dot(v2); Real dot11 = v1.dot(v1); Real dot12 = v1.dot(v2); // Compute barycentric coordinates Real invDenom = 1.0 / (dot00 * dot11 - dot01 * dot01); u = (dot11 * dot02 - dot01 * dot12) * invDenom; v = (dot00 * dot12 - dot01 * dot02) * invDenom; // Check if point is in triangle return (u >= -1e-8) && (v >= -1e-8) && (u + v < 1.0 + 1e-8); } /* -------------------------------------------------------------------------- */ // returns the number of the element an atom is in // or -1 if it's outside the template inline int ComputeCADDTemplate::findElement(Vector &at_coords) { auto &nodes = this->mesh.getContainerNodes(); int i = 0; for (auto &&elem : this->mesh.getContainerElems()) { auto node1 = nodes[elem.localIndexes()[0]].position0(); auto node2 = nodes[elem.localIndexes()[1]].position0(); auto node3 = nodes[elem.localIndexes()[2]].position0(); if (is_in_triangle(node1, node2, node3, at_coords.block<2, 1>(0, 0))) { return i; } ++i; } return UINT_MAX; } /* -------------------------------------------------------------------------- */ void ComputeCADDTemplate::init() { loadNPZ(filename); } /* -------------------------------------------------------------------------- */ /* LMDESC CADD_TEMPLATE TODO */ /* LMEXAMPLE COMPUTE template CADD_TEMPLATE INPUT md TODO */ void ComputeCADDTemplate::declareParams() { ComputeInterface::declareParams(); /* LMKEYWORD FILENAME Set the filename containing a core template */ this->parseKeyword("FILENAME", filename); // /* LMKEYWORD ANGLE_BETWEEN_B_AND_X // Specifies the angle between Burgers vector and X(11-2) axis // */ // this->parseKeyword("ANGLE_BETWEEN_B_AND_X", angle_between_b_and_x); } /* -------------------------------------------------------------------------- */ -void ComputeCADDTemplate::coreDispAtom(Vector &pos_core, - Vector &disp) { +template +void ComputeCADDTemplate::coreDispAtom(Pos &&pos_core, Disp &&disp) { int elem_num = this->findElement(pos_core); + if (elem_num == -1) + return; - if (elem_num != -1) { - Vector<3> coords = {1.0, pos_core[0], pos_core[1]}; - auto &elems = this->mesh.getContainerElems(); - auto &nodes = this->mesh.getContainerNodes(); - auto &elem = elems[elem_num]; - - auto &node1 = nodes[elem.localIndexes()[0]].position0(); - auto &node2 = nodes[elem.localIndexes()[1]].position0(); - auto &node3 = nodes[elem.localIndexes()[2]].position0(); - - Eigen::Matrix Su; - Su << 1., 0., 0., node1, node2, node3; - Matrix<3> SxInv; - auto &&tmp = Su * SxInv * coords; - disp += tmp.block<3, 1>(1, 0); - } + Vector<3> coords = {1.0, pos_core[0], pos_core[1]}; + auto &elems = this->mesh.getContainerElems(); + auto &nodes = this->mesh.getContainerNodes(); + auto &elem = elems[elem_num]; + + auto &node1 = nodes[elem.localIndexes()[0]].position0(); + auto &node2 = nodes[elem.localIndexes()[1]].position0(); + auto &node3 = nodes[elem.localIndexes()[2]].position0(); + + auto &u1 = this->corrections[elem.localIndexes()[0]]; + auto &u2 = this->corrections[elem.localIndexes()[1]]; + auto &u3 = this->corrections[elem.localIndexes()[2]]; + + Eigen::Matrix Sx = Eigen::Matrix::Zero(); + + Sx.block<1, 3>(0, 0) = Vector<3>{1., 0., 0.}; + Sx.block<2, 1>(1, 0) = node3; + Sx.block<2, 1>(1, 1) = node1 - node3; + Sx.block<2, 1>(1, 2) = node2 - node3; + + Matrix<3> SxInv = Sx.inverse(); + + Eigen::Matrix Su = Eigen::Matrix::Zero(); + + Su.block<1, 3>(0, 0) = Vector<3>{1., 0., 0.}; + Su.block<3, 1>(1, 0) = u3; + Su.block<3, 1>(1, 1) = u1 - u3; + Su.block<3, 1>(1, 2) = u2 - u3; + + auto &&tmp = Su * SxInv * coords; + Vector<3> _disp = tmp.block<3, 1>(1, 0); + disp += _disp; } /* -------------------------------------------------------------------------- */ template std::enable_if_t ComputeCADDTemplate::build(Points &points, Segment &segment) { - for (auto &&point : points) { - computeTemplateDisplacements(point, segment); + auto &displacements = this->getOutputAsArray(); + displacements.resize(points.size(), Dim); + + for (auto &&[point, disp] : zip(points, displacements.rowwise())) { + Vector<3> _disp; + computeTemplateDisplacements(point, segment, _disp); + disp = _disp; } } /* -------------------------------------------------------------------------- */ template Matrix<3u> computeTemplateRotation(Vec1 &&_dd_line_dir, Vec2 &&_slip_normal_dir) { Vector<3u> dd_line_dir(_dd_line_dir); Vector<3u> slip_normal_dir(_slip_normal_dir); // construct a rotation where: // - the normal to slip plane becomes z // - the line direction becomes y // - normal to the line in the slip plane becomes x - // checks is the normal is normal to the line direction + // checks if computed normal is orthogonal to the line direction LM_ASSERT(std::abs(dd_line_dir.dot(slip_normal_dir)) < 1e-10, "normal and line not defined properly"); // construct the first basis vector Vector<3u> first_basis_vector = dd_line_dir.cross(slip_normal_dir); Matrix<3u> Rmat = Matrix<3>::Zero(); Rmat.block<3, 1>(0, 0) = first_basis_vector; Rmat.block<3, 1>(0, 1) = dd_line_dir; Rmat.block<3, 1>(0, 2) = slip_normal_dir; return Rmat.transpose(); } /* -------------------------------------------------------------------------- */ -template -void ComputeCADDTemplate::computeTemplateDisplacements(Point &p, - ElementType &segments) { +template +void ComputeCADDTemplate::computeTemplateDisplacements(Point &&p, + ElementType &&segments, + Disp &&disp) { for (auto &&segment : segments.getContainerElems()) { auto &&burgers = segment.burgers(); auto &&normal = segment.normal(); LM_ASSERT(std::abs(burgers.norm()) > 1e-8, "No found burgers vector for neighbors"); auto &&pos0 = p.position0(); auto conn = segment.localIndexes(); auto node1 = segment.getNode(0).position(); auto node2 = segment.getNode(1).position(); // computes the direction of the line auto dd_line = node2 - node1; // computes a unitary vector in the direction of the line Vector dd_line_direction = dd_line; dd_line_direction.normalize(); // computes the relatice position of the atom Vector pos_atom = pos0 - node1; // Rmat is a rotation of only the coordinates in the slip plane coordinates // This is to be used to convert to/from the reference where the templates // were computed Matrix Rmat = computeTemplateRotation(dd_line_direction, normal); // project the atom position to the segment Vector projected_atom_pos = pos_atom.dot(dd_line_direction) * dd_line_direction; // vector from segment to atom position Vector vec_atom_2_dd_core = pos_atom - projected_atom_pos; Vector template_coords = Rmat * vec_atom_2_dd_core; - Vector disp = Vector::Zero(); this->coreDispAtom(template_coords, disp); - LM_TOIMPLEMENT; - // for (UInt i = 0; i < Dim; ++i) - // this->push_back(disp[i]); } } /* -------------------------------------------------------------------------- */ //**************** // part concrening the interpolation between angles //**************** // character angle: angle between burgers vector and dislocation line // Real dangle = fabs(burgers.dot(dd_line_direction)) / burgers.norm(); // Real dangle_deg = // std::acos(dangle) * 180.0 / M_PI; // character angle (degree) // Real distance = fabs(90.0 - dangle_deg); // Real cangle = distance + 90.0; // change in a range between 90 and 180 // Real closest_angles[2] = {0.0, 0.0}; // UInt indices_closest_angles[2] = {1, 1}; // compute_angle_ratio(cangle, indices_closest_angles, closest_angles); // DO NOT DELETE: displace elsewhere !!!!! // contained the displacement of the two templates // involved for angle interpolation // std::vector> d = {0.0, 0.0}; // DO NOT DELETE: displace elsewhere !!!!! // loop two neighboring angle templates for interpolation // for (int j = 0; j < 2; ++j) { // Real tmp_angle = fabs(this->angle_between_b_and_x - cangle); // tdisp = Rmat.transpose() * disp_template0; // tdisp[0] = disp_template0[2] * cos(tmp_angle / 180.0 * M_PI) + // disp_template0[0] * sin(tmp_angle / 180.0 * M_PI); // tdisp[1] = disp_template0[1]; // tdisp[2] = disp_template0[0] * cos(tmp_angle / 180.0 * M_PI) - // disp_template0[2] * sin(tmp_angle / 180.0 * M_PI); // DO NOT DELETE: displace elsewhere !!!!! // in between angles interpolation // tdisp[0] = (dx[1] - dx[0]) / (nei_angles[1] - nei_angles[0]) * // (cangle - nei_angles[0]) + // dx[0]; // tdisp[1] = (dy[1] - dy[0]) / (nei_angles[1] - nei_angles[0]) * // (cangle - nei_angles[0]) + // dy[0]; // tdisp[2] = (dz[1] - dz[0]) / (nei_angles[1] - nei_angles[0]) * // (cangle - nei_angles[0]) + // dz[0]; // /* -------------------------------------------------------------------------- // */ // template // bool CADD_template::PointInSquare( // Vector &pts1, Vector &pts2, Vector &pts3, Vector // &pts4, // Vector &pts, UInt slip_case) { // Real area0 = 0.0; // Real area1 = AreaTriangle(pts1, pts3, pts4, slip_case); // Real area2 = AreaTriangle(pts4, pts2, pts1, slip_case); // Real area3 = AreaTriangle(pts1, pts2, pts3, slip_case); // Real area4 = AreaTriangle(pts2, pts3, pts4, slip_case); // area0 = (fabs(area1) + fabs(area2) + fabs(area3) + fabs(area4)) / 2.; // Vector<2> int_X; // if (intersectTwoD(pts1, pts3, pts2, pts4, int_X, slip_case)) { // Real temp[3] = {0.0, 0.0, 0.0}; // for (UInt i = 0; i < 3; ++i) { // temp[i] = pts2[i]; // pts2[i] = pts1[i]; // pts1[i] = temp[i]; // } // } // Real tot_area = 0.0; // Real tot_area1 = AreaTriangle(pts, pts1, pts2, slip_case); // Real tot_area2 = AreaTriangle(pts, pts3, pts1, slip_case); // Real tot_area3 = AreaTriangle(pts, pts4, pts3, slip_case); // Real tot_area4 = AreaTriangle(pts, pts2, pts4, slip_case); // tot_area = // fabs(tot_area1) + fabs(tot_area2) + fabs(tot_area3) + fabs(tot_area4); // if (area0 + 1.0 < tot_area) // return false; // else // return true; // } // template void ComputeCADDTemplate<_Input>::init() { // if (this->is_in_atomic) { // for (UInt i = 0; i < this->num_templates; ++i) { // // TwoDTemplate *single_template = // // new TwoDTemplate(); // // this->core_templates.push_back(single_template); // // this->core_templates[i]->init(this->templatefilename_list[i]); // // } // LM_TOIMPLEMENT; // } // } // this function retreives the two closest angles // this is done to make the interpolation between angles easier // TODO: review move elsewhere inline void compute_angle_ratio(Real angle [[gnu::unused]], UInt *indices_closest_angles [[gnu::unused]], Real *closest_angles [[gnu::unused]]) { // TO REVIEW; // to enforce the specific core template // if(this->use_specific_template) { // angle = this->specific_template_angle; // } else { // if (angle < 90.0) { // angle = 180.0 - angle; // } else if ((180.0 < angle) && (angle <= 270.0)) { // angle = 360.0 - angle; // } else if ((270.0 < angle) && (angle <= 360.0)) { // angle = angle - 180.0; // } // // } // for (UInt j = 0; j < this->num_templates - 1; ++j) { // two_angles[0] = this->core_templates[j]->get_theta(); // two_angles[1] = this->core_templates[j + 1]->get_theta(); // two_indices[0] = j; // two_indices[1] = j + 1; // if ((angle - two_angles[0]) * (angle - two_angles[1]) <= 0.0) { // break; // } // } } /* -------------------------------------------------------------------------- */ void ComputeCADDTemplate::compute_make_call() { DUMP(this->getID() << ": Compute", DBG_INFO); this->build(this->getInput("points"), this->getInput("segments")); } /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ diff --git a/src/filter/compute_cadd_template.hh b/src/filter/compute_cadd_template.hh index 3a7653f..f3350f4 100644 --- a/src/filter/compute_cadd_template.hh +++ b/src/filter/compute_cadd_template.hh @@ -1,121 +1,123 @@ /** * @file two2_template.hh * * @author Guillaume Anciaux * @author Till Junge * @author Jaehyun Cho * * @date Wed Jan 15 17:00:43 2014 * * @brief 2d template helper * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * LibMultiScale is free software: you can redistribute it and/or modify it * under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * LibMultiScale is distributed in the hope that it will be useful, but * WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with LibMultiScale. If not, see . * */ #ifndef __LIBMULTISCALE_COMPUTE_CADD_TEMPLATE_HH__ #define __LIBMULTISCALE_COMPUTE_CADD_TEMPLATE_HH__ /* ------------------------------------------------------------------------ */ #include "compute_interface.hh" #include "container_mesh.hh" #include "ref_elem_dd.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ /* the use_template switch turns on or off the use of the template for the dislocation detection, not for the coupling */ class ComputeCADDTemplate : public ComputeInterface { public: static constexpr UInt Dim = 3; DECLARE_COMPUTE(ComputeCADDTemplate, _or, _or); template std::enable_if_t build(Points &points, Segment &segments); template std::enable_if_t build(Points &, Segment &) { LM_TOIMPLEMENT; } void init() override; private: void loadNPZ(std::string &filename); - template - void computeTemplateDisplacements(Point &point, ElementType &segment); + template + void computeTemplateDisplacements(Point &&point, ElementType &&segment, + Disp &&disp); - void coreDispAtom(Vector &pos_core, Vector &disp); + template + void coreDispAtom(Pos &&pos_core, Disp &&disp); template Real evaluateSquareError(const Real pos_core[], IteratorAtoms &it_atoms); Real get_r_cut() { return this->r_cut; } Real get_burgers() { return this->burgers; } Real get_theta() { return this->theta; } inline UInt get_search_dim() { return this->search_dim; } inline Real get_tol() { return this->tol; } Real domainSize[3]; private: inline int findElement(Vector &at_coords); ContainerGenericMesh<2> mesh; std::vector> corrections; Real center_atoms_around_core[2]; Real r_cut; UInt num_nodes, num_elems; std::vector> shape; // conversion between xy coordinates and element // coordinates l1, l2, l3 std::vector> x_l_transition; Vector md_bbox_min, md_bbox_max; Real poisson, burgers, theta; dislo_orientation inversion; Real factor[2]; // depends entirely on inversion // some temp data // Real atom_coords[3]; //! Optimization parameters and temps UInt search_dim; Real tol; // Real **vertices, *vertices_cont; // UInt maxiter; // int dir_replica; // int nb_replica; // Real rcut; std::string filename; Real angle_between_b_and_x; }; /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ #endif diff --git a/src/filter/compute_dd_mesh.cc b/src/filter/compute_dd_mesh.cc index 3b9b542..0bb9df5 100644 --- a/src/filter/compute_dd_mesh.cc +++ b/src/filter/compute_dd_mesh.cc @@ -1,128 +1,129 @@ /** * @file compute_mesh.cc * * @author Guillaume Anciaux * * @date Wed Jul 09 21:59:47 2014 * * @brief This computes generate an independent mesh from points and * connectivities provided as computes * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * LibMultiScale is free software: you can redistribute it and/or modify it * under the * terms of the GNU Lesser General Public License as published by the * Free * Software Foundation, either version 3 of the License, or (at your option) * any * later version. * * LibMultiScale is distributed in the hope that it will be useful, but * WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for * more * details. * * You should have received a copy of the GNU Lesser General Public * License * along with LibMultiScale. If not, see . * */ /* -------------------------------------------------------------------------- */ // #include "compute_interface.hh" // #include "compute_point_set.hh" // #include "container_mesh.hh" #include "factory_multiscale.hh" #include "lib_continuum.hh" #include "lib_dd.hh" #include "lib_dumper.hh" #include "lib_filter.hh" #include "lib_md.hh" #include "lib_stimulation.hh" #include "lm_common.hh" #include "lm_communicator.hh" /* ----------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* ----------------------------------------------------------------------- */ ComputeDDMesh::ComputeDDMesh(const std::string &name) : LMObject(name), FilterInterface() { this->createInput("nodes"); this->createInput("connectivity"); this->createInput("burgers"); this->createInput("normals"); this->createOutput("mesh"); } /* ----------------------------------------------------------------------- */ ComputeDDMesh::~ComputeDDMesh() {} /* ----------------------------------------------------------------------- */ template void ComputeDDMesh::build(Nodes &nodes, Connectivity &connectivity, Burgers &burgers, Normals &normals) { constexpr UInt Dim = Nodes::Dim; using MeshType = ContainerGenericDDMesh; auto &output = allocOutput("mesh", this->getID()); auto &nds = output.getContainerNodes(); nds.clear(); for (auto &&nd : nodes) { RefGenericDDNode ref = nd; nds.push_back(ref); } auto &els = output.getContainerElems(); els.clear(); - auto nodes_per_elem = connectivity.getDim(); - auto size = connectivity.size(); - auto nb_elem = size / nodes_per_elem; RefGenericDDElem tmp_el; - for (UInt i = 0; i < nb_elem; ++i) { - Eigen::VectorXd conn(nodes_per_elem); - for (UInt j = 0; j < nodes_per_elem; ++j) { - conn[j] = UInt(connectivity[i * nodes_per_elem + j]); - tmp_el.setConnectivity(conn); - } + + for (auto &&[node_indexes, burg, norm] : + zip(connectivity.rowwise(), burgers.rowwise(), normals.rowwise())) { + tmp_el.setConnectivity(node_indexes); + tmp_el.setNodes(nds); + tmp_el.burgers() = burg; + tmp_el.normal() = norm; els.push_back(tmp_el); } } -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ /* LMDESC MESH This compute generates an independent mesh from an input containing points and a compute which contains connectivities. */ -/* LMEXAMPLE COMPUTE pset DDMESH INPUT nodes=coords INPUT connectivity=conn +/* LMEXAMPLE COMPUTE pset DDMESH INPUT nodes=coords INPUT connectivity=conn INPUT burgers=burg INPUT normals=ns */ void ComputeDDMesh::declareParams() { FilterInterface::declareParams(); } -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ void ComputeDDMesh::compute_make_call() { this->build(this->getInput("nodes"), this->getInput("connectivity"), this->getInput("burgers"), this->getInput("normals")); } -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ __END_LIBMULTISCALE__ diff --git a/src/filter/compute_impulse.cc b/src/filter/compute_impulse.cc index 02879ac..4317ddc 100644 --- a/src/filter/compute_impulse.cc +++ b/src/filter/compute_impulse.cc @@ -1,202 +1,202 @@ /** * @file compute_impulse.cc * * @author Guillaume Anciaux * * @date Wed Jul 09 21:59:47 2014 * * @brief This stimulator sets an initial impulse displacement field that will * generate a wave * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * LibMultiScale is free software: you can redistribute it and/or modify it * under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * LibMultiScale is distributed in the hope that it will be useful, but * WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with LibMultiScale. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "compute_impulse.hh" #include "lib_continuum.hh" #include "lib_dd.hh" #include "lib_md.hh" #include "lm_common.hh" #include "ref_point_data.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ ComputeImpulse::ComputeImpulse(const std::string &name) : LMObject(name) { impulseType = Impulse::PLANAR; direction = 0; center.setZero(); this->createInput("domain"); this->createOutput("output") = ContainerArray(this->getID() + ":" + "output"); } /* -------------------------------------------------------------------------- */ ComputeImpulse::~ComputeImpulse() {} /* -------------------------------------------------------------------------- */ void ComputeImpulse::init() { wave_number = 2 * M_PI / wave_length; } /* -------------------------------------------------------------------------- */ inline Real ComputeImpulse::computeScaleImpulse(Real d) { const Real k = wave_number; const Real res = intensity * sin(k * d) * exp(-k * k / 9 * d * d); DUMP("impulse computed: " << d << " " << k << " " << intensity << " => " << res, DBG_INFO); return res; } /* -------------------------------------------------------------------------- */ template <> inline UInt Parser::parse(Impulse::ImpulseType &type, std::stringstream &line, UInt) { std::string name; UInt nb = strNext(name, line); if (name == "CIRCULAR") type = Impulse::CIRCULAR; else if (name == "PLANAR") type = Impulse::PLANAR; else LM_FATAL("unknown impulse type " << name); return nb; } /* -------------------------------------------------------------------------- */ template void ComputeImpulse::build(Cont &cont) { constexpr UInt Dim = Cont::Dim; this->clear(); auto &output = this->getOutputAsArray(); + output.resize(0, Dim); for (auto &&at : cont) { if (impulseType == Impulse::PLANAR) { Real d = at.position0()[direction] - center[direction]; Real impulse = computeScaleImpulse(d); for (UInt i = 0; i < Dim; ++i) { if (i == direction) output.push_back(impulse); } } else if (impulseType == Impulse::CIRCULAR) { - output.resize(output.rows(), Dim); Vector vec_dir = at.position0() - center.cast().block(0, 0); Real norm_vec_dir = vec_dir.norm(); if (norm_vec_dir == 0.) LM_FATAL("Cannot compute a displacement at the center of the circle"); vec_dir /= norm_vec_dir; Real d = norm_vec_dir - radius; Real impulse = computeScaleImpulse(d); - this->getOutputAsArray().push_back(impulse * vec_dir); + output.push_back(impulse * vec_dir); } } } /* -------------------------------------------------------------------------- */ /* LMDESC IMPULSE This stimulator sets an initial impulse displacement field that will generate a wave. It can generate planar or circular waves: - Planar: .. math:: u(\vec{X}) += I \cdot \sin\left(k d \right) \exp\left(-\frac{(k d)^2}{9} \right) \vec{e_i} with :math:`I` the intensity, :math:`i` the chosen direction and :math:`k = 2\pi/LWAVE` and :math:`d = \vec{X}_i-\vec{C}_i` as defined with the parameters (see below). - Circular: .. math:: u(X) += I \cdot \sin\left(k (d - R) \right) \exp\left(-\frac{(k (d - R))^2}{9} \right) \frac{\vec{X} - \vec{C}}{d} with :math:`I` the intensity, :math:`k = 2\pi/LWAVE` and :math:`d = \mid \vec{X} - \vec{C} \mid` and :math:`R` a radius as defined with the parameters (see below). */ /* LMEXAMPLE COMPUTE impulse IMPULSE INPUT central_md LWAVE WaveLength*r0 * DIRECTION 0 INTENSITY 1.1e-3 ONESHOT 0 */ void ComputeImpulse::declareParams() { /* LMKEYWORD LWAVE Set the wavelength of the wave to be generated */ this->parseKeyword("LWAVE", wave_length); /* LMKEYWORD DIRECTION It can be 0,1 or 2 to specify X, Y or Z direction for planar waves */ this->parseKeyword("DIRECTION", direction, 0u); /* LMKEYWORD RADIUS Radius to be used to generate the impulse */ this->parseKeyword("RADIUS", radius, 0.); /* LMKEYWORD INTENSITY Intensity of the displacement */ this->parseKeyword("INTENSITY", intensity); /* LMKEYWORD CENTER center of the impulse to be injected */ this->parseVectorKeyword("CENTER", spatial_dimension, center, VEC_DEFAULTS(0., 0., 0.)); /* LMKEYWORD TYPE Flag to specify the type of impulse. It can be either CIRCULAR or PLANAR. Self-explanatory examples are: - COMPUTE impulse-planar IMPULSE INPUT md LWAVE L0/10 DIRECTION 0 INTENSITY 1e-3 CENTER 0 0 .. image:: images/planar-impulse.png - COMPUTE impulse-circular IMPULSE INPUT md LWAVE L0/10 TYPE CIRCULAR RADIUS L0/4 INTENSITY 1e-3 CENTER 0 0 .. image:: images/radial-impulse.png */ this->parseKeyword("TYPE", impulseType, Impulse::PLANAR); } /* -------------------------------------------------------------------------- */ DECLARE_COMPUTE_MAKE_CALL(ComputeImpulse) __END_LIBMULTISCALE__