Author: andrico Date: Fri Nov 4 17:32:17 2016 New Revision: 5463 URL: http://svn.gna.org/viewcvs/getfem?rev=5463&view=rev Log: interpolate_transformation_on_deformed_domains. work in progress
Added: trunk/getfem/src/getfem_interpolation_on_deformed_domains.cpp (with props) Modified: trunk/getfem/msvc/libgetfem/libgetfem.vcxproj trunk/getfem/src/getfem/getfem_generic_assembly.h Modified: trunk/getfem/msvc/libgetfem/libgetfem.vcxproj URL: http://svn.gna.org/viewcvs/getfem/trunk/getfem/msvc/libgetfem/libgetfem.vcxproj?rev=5463&r1=5462&r2=5463&view=diff ============================================================================== --- trunk/getfem/msvc/libgetfem/libgetfem.vcxproj (original) +++ trunk/getfem/msvc/libgetfem/libgetfem.vcxproj Fri Nov 4 17:32:17 2016 @@ -271,6 +271,7 @@ <ClCompile Include="..\..\src\getfem_integration_composite.cc" /> <ClCompile Include="..\..\src\getfem_interpolated_fem.cc" /> <ClCompile Include="..\..\src\getfem_interpolation.cc" /> + <ClCompile Include="..\..\src\getfem_interpolation_on_deformed_domains.cpp" /> <ClCompile Include="..\..\src\getfem_level_set.cc" /> <ClCompile Include="..\..\src\getfem_level_set_contact.cc" /> <ClCompile Include="..\..\src\getfem_linearized_plates.cc" /> Modified: trunk/getfem/src/getfem/getfem_generic_assembly.h URL: http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_generic_assembly.h?rev=5463&r1=5462&r2=5463&view=diff ============================================================================== --- trunk/getfem/src/getfem/getfem_generic_assembly.h (original) +++ trunk/getfem/src/getfem/getfem_generic_assembly.h Fri Nov 4 17:32:17 2016 @@ -552,6 +552,27 @@ (ga_workspace &workspace, const std::string &transname, const mesh &source_mesh, const mesh &target_mesh, const std::string &expr); + /** Add a transformation to the workspace that creates an identity mapping between + two meshes in deformed state. Conceptually, it can be viewed as a transformation + from expression Xsource + Usource - Utarget, except such an expression + cannot be used directly in the transformation from expression (function above), + as Utarget needs to be interpolated though an inversion of the transformation of + the target domain. + Thread safe if added to thread local workspace. + */ + void add_interpolate_transformation_on_deformed_domains + (ga_workspace &workspace, const std::string &transname, + const mesh &source_mesh, const std::string &source_displacements, + const mesh_region &source_region, const mesh &target_mesh, + const std::string &target_displacements, const mesh_region &target_region); + + /**.. the same as above, but adding transformation to the model. + Note, this version is not thread safe.*/ + void add_interpolate_transformation_on_deformed_domains + (model &md, const std::string &transname, + const mesh &source_mesh, const std::string &source_displacements, + const mesh_region &source_region, const mesh &target_mesh, + const std::string &target_displacements, const mesh_region &target_region); /** Create a new instance of a transformation corresponding to the interpolation on the neighbour element. Can only be applied to the Added: trunk/getfem/src/getfem_interpolation_on_deformed_domains.cpp URL: http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_interpolation_on_deformed_domains.cpp?rev=5463&view=auto ============================================================================== --- trunk/getfem/src/getfem_interpolation_on_deformed_domains.cpp (added) +++ trunk/getfem/src/getfem_interpolation_on_deformed_domains.cpp Fri Nov 4 17:32:17 2016 @@ -0,0 +1,362 @@ +/*=========================================================================== + + Copyright (C) 2013-2016 Yves Renard, Konstantinos Poulios and Andriy Andreykiv. + + This file is a part of GetFEM++ + + GetFEM++ 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 along with the GCC Runtime Library + Exception either version 3.1 or (at your option) any later version. + This program 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 and GCC Runtime Library Exception for more details. + You should have received a copy of the GNU Lesser General Public License + along with this program; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. + +===========================================================================*/ + +#include "getfem/getfem_generic_assembly.h" + +namespace getfem { + +// Structure describing a contact boundary (or contact body) +struct contact_boundary { + size_type region; // boundary region for the slave (source) + // and volume region for the master (target) + const getfem::mesh_fem *mfu; // F.e.m. for the displacement. + std::string dispname; // Variable name for the displacement + mutable const model_real_plain_vector *U;// Displacement + mutable model_real_plain_vector U_unred; // Unreduced displacement + + contact_boundary(size_type r, const mesh_fem *mf, const std::string &dn) + : region(r), mfu(mf), dispname(dn) {} +}; + +base_small_vector element_U(const contact_boundary &cb, size_type cv) +{ + auto U_elm = base_small_vector{}; + slice_vector_on_basic_dof_of_element(*(cb.mfu), *cb.U, cv, U_elm); + return U_elm; +} + +//Returns an iterator of a box which centre is closest to the given point +auto most_central_box(const bgeot::rtree::pbox_set &bset, + const bgeot::base_node &pt) -> decltype(begin(bset)) +{ + using namespace std; + + auto itmax = begin(bset); + + auto it = itmax; + if (bset.size() > 1) { + auto rate_max = scalar_type{-1}; + for (; it != bset.end(); ++it) { + auto rate_box = scalar_type{1}; + for (size_type i = 0; i < pt.size(); ++i) { + auto h = (*it)->max[i] - (*it)->min[i]; + if (h > 0.) { + auto rate = min((*it)->max[i] - pt[i], pt[i] - (*it)->min[i]) / h; + rate_box = min(rate, rate_box); + } + } + if (rate_box > rate_max) { + itmax = it; + rate_max = rate_box; + } + } + } + + return itmax; +} + +class interpolate_transformation_on_deformed_domains + : public virtual_interpolate_transformation { + + contact_boundary master; + contact_boundary slave; + + mutable bgeot::rtree element_boxes; + mutable std::vector<size_type> box_to_convex; //index to obtain + //a convex number from a box number + mutable bgeot::geotrans_inv_convex gic; + mutable fem_precomp_pool fppool; + + //Create a box tree based on the deformed elements of the master (target) + void compute_element_boxes() const { // called by init + base_matrix G; + model_real_plain_vector Uelm; //element displacement + element_boxes.clear(); + + auto bnum = master.region; + auto &mfu = *(master.mfu); + auto &U = *(master.U); + auto &m = mfu.linked_mesh(); + auto N = m.dim(); + + base_node Xdeformed(N), bmin(N), bmax(N); + auto region = m.region(bnum); + + //the box tree creation and subsequent transformation inversion + //should be done for all elements of the master, while integration + //will be performed only on a thread partition of the slave + region.prohibit_partitioning(); + + GMM_ASSERT1(mfu.get_qdim() == N, "Wrong mesh_fem qdim"); + + dal::bit_vector points_already_interpolated; + std::vector<base_node> transformed_points(m.nb_max_points()); + box_to_convex.clear(); + box_to_convex.reserve(region.size()); + + for (getfem::mr_visitor v(region, m); !v.finished(); ++v) { + auto cv = v.cv(); + auto pgt = m.trans_of_convex(cv); + auto pf_s = mfu.fem_of_element(cv); + auto pfp = fppool(pf_s, pgt->pgeometric_nodes()); + + slice_vector_on_basic_dof_of_element(mfu, U, cv, Uelm); + mfu.linked_mesh().points_of_convex(cv, G); + + auto ctx = fem_interpolation_context{pgt, pfp, size_type(-1), G, cv}; + auto nb_pt = pgt->structure()->nb_points(); + + for (size_type k = 0; k < nb_pt; ++k) { + auto ind = m.ind_points_of_convex(cv)[k]; + + // computation of a transformed vertex + ctx.set_ii(k); + if (points_already_interpolated.is_in(ind)) { + Xdeformed = transformed_points[ind]; + } else { + pf_s->interpolation(ctx, Uelm, Xdeformed, dim_type{N}); + Xdeformed += ctx.xreal(); //Xdeformed = U + Xref + transformed_points[ind] = Xdeformed; + points_already_interpolated.add(ind); + } + + if (k == 0) // computation of bounding box + bmin = bmax = Xdeformed; + else { + for (size_type l = 0; l < N; ++l) { + bmin[l] = std::min(bmin[l], Xdeformed[l]); + bmax[l] = std::max(bmax[l], Xdeformed[l]); + } + } + } + + // Safety coefficient of 1.3 (for nonlinear transformations) + // Expanding the box by 15% of the largest edge in every direction + auto h = bmax[0] - bmin[0]; + for (size_type k = 1; k < N; ++k) h = std::max(h, bmax[k] - bmin[k]); + for (size_type k = 0; k < N; ++k) {bmin[k] -= h * 0.15; bmax[k] += h * 0.15;} + + // Store the bounding box and additional information. + element_boxes.add_box(bmin, bmax, box_to_convex.size()); + box_to_convex.push_back(cv); + } + } + +public: + + interpolate_transformation_on_deformed_domains( + size_type source_region, + const getfem::mesh_fem &mf_source, + const std::string &source_displacements, + size_type target_region, + const getfem::mesh_fem &mf_target, + const std::string &target_displacements) + : + slave{source_region, &mf_source, source_displacements}, + master{target_region, &mf_target, target_displacements} +{} + + + void extract_variables(const ga_workspace &workspace, + std::set<var_trans_pair> &vars, + bool ignore_data, + const mesh &m_x, + const std::string &interpolate_name) const override { + if (!ignore_data || !(workspace.is_constant(master.dispname))){ + vars.emplace(master.dispname, interpolate_name); + vars.emplace(slave.dispname, ""); + } + } + + void init(const ga_workspace &workspace) const override { + + for (auto pcb : {&master, &slave}) { + auto &mfu = *(pcb->mfu); + if (mfu.is_reduced()) { + gmm::resize(pcb->U_unred, mfu.nb_basic_dof()); + mfu.extend_vector(workspace.value(pcb->dispname), pcb->U_unred); + pcb->U = &(pcb->U_unred); + } else { + pcb->U = &(workspace.value(pcb->dispname)); + } + } + + compute_element_boxes(); + }; + + void finalize() const override { + element_boxes.clear(); + box_to_convex.clear(); + master.U_unred.clear(); + slave.U_unred.clear(); + fppool.clear(); + } + + std::vector<bgeot::base_node> deformed_master_nodes(size_type cv) const { + using namespace bgeot; + using namespace std; + + auto nodes = vector<base_node>{}; + + auto U_elm = element_U(master, cv); + auto &mfu = *(master.mfu); + auto G = base_matrix{}; + auto pfu = mfu.fem_of_element(cv); + auto pgt = master.mfu->linked_mesh().trans_of_convex(cv); + auto pfp = fppool(pfu, pgt->pgeometric_nodes()); + auto N = mfu.linked_mesh().dim(); + auto pt = base_node(N); + auto U = base_small_vector(N); + master.mfu->linked_mesh().points_of_convex(cv, G); + auto ctx = fem_interpolation_context{pgt, pfp, size_type(-1), G, cv}; + auto nb_pt = pgt->structure()->nb_points(); + nodes.reserve(nb_pt); + for (size_type k = 0; k < nb_pt; ++k) { + ctx.set_ii(k); + pfu->interpolation(ctx, U_elm, U, dim_type{N}); + gmm::add(ctx.xreal(), U, pt); + nodes.push_back(pt); + } + + return nodes; + } + + int transform(const ga_workspace &workspace, + const mesh &m_x, + fem_interpolation_context &ctx_x, + const base_small_vector &/*Normal*/, + const mesh **m_t, + size_type &cv, + short_type &face_num, + base_node &P_ref, + base_small_vector &N_y, + std::map<var_trans_pair, base_tensor> &derivatives, + bool compute_derivatives) const override { + + auto &target_mesh = master.mfu->linked_mesh(); + *m_t = &target_mesh; + int ret_type = 0; + + using namespace gmm; + using namespace bgeot; + using namespace std; + + //compute a deformed point of the slave (also using terminology: source or x) + auto cv_x = ctx_x.convex_num(); + auto U_elm_x = element_U(slave, cv_x); + auto &mfu_x = *(slave.mfu); + auto pfu_x = mfu_x.fem_of_element(cv_x); + auto N = mfu_x.linked_mesh().dim(); + auto U_x = base_small_vector(N); + auto G_x = base_matrix{}; //coordinates of the source element nodes + m_x.points_of_convex(cv_x, G_x); + ctx_x.set_pf(pfu_x); + pfu_x->interpolation(ctx_x, U_elm_x, U_x, dim_type{N}); + auto pt_x = base_small_vector(N); //deformed point of the slave + add(ctx_x.xreal(), U_x, pt_x); + + //Find the best box from the master (target) that + //corresponds to this point (The box which centre is the closest to the point). + //Obtain the corresponding element number using the box id and box_to_convex + //indices. Compute deformed nodes of the target element. Invert the geometric + //transformation of the target element with deformed nodes, obtaining this way + //reference coordinates of the target element + auto bset = rtree::pbox_set{}; + element_boxes.find_boxes_at_point(pt_x, bset); + while (!bset.empty()) + { + auto itmax = most_central_box(bset, pt_x); + cv = box_to_convex[(*itmax)->id]; + auto deformed_nodes_y = deformed_master_nodes(cv); + gic.init(deformed_nodes_y, target_mesh.trans_of_convex(cv)); + auto converged = true; + auto is_in = gic.invert(pt_x, P_ref, converged, 1E-4); + if (is_in && converged) { + face_num = static_cast<short_type>(-1); + ret_type = 1; + break; + } + if (bset.size() == 1) break; + bset.erase(itmax); + } + + //Since this transformation can be seen as Xsource + Usource - Utarget, + //the corresponding stiffnesses are identity matrix for Usource and + //minus identity for Utarget + if (compute_derivatives && ret_type == 1) { + GMM_ASSERT2(derivatives.size() == 2, + "Expecting to return derivatives only for Umaster and Uslave"); + for (auto &pair : derivatives) + { + if (pair.first.varname == slave.dispname) + { + for (size_type i = 0; i < m_x.dim(); ++i) pair.second(i, i) = 1.0; + } + else + if (pair.first.varname == master.dispname) + { + for (size_type i = 0; i < target_mesh.dim(); ++i) pair.second(i, i) = -1.0; + } + else GMM_ASSERT2(false, "unexpected derivative variable"); + } + } + return ret_type; + } + +}; + + void add_interpolate_transformation_on_deformed_domains + (ga_workspace &workspace, const std::string &transname, + const mesh &source_mesh, const std::string &source_displacements, + const mesh_region &source_region, const mesh &target_mesh, + const std::string &target_displacements, const mesh_region &target_region) + { + auto pmf_source = workspace.associated_mf(source_displacements); + auto pmf_target = workspace.associated_mf(target_displacements); + auto p_transformation + = std::make_shared<interpolate_transformation_on_deformed_domains>(source_region.id(), + *pmf_source, + source_displacements, + target_region.id(), + *pmf_target, + target_displacements); + workspace.add_interpolate_transformation(transname, p_transformation); + } + + void add_interpolate_transformation_on_deformed_domains + (model &md, const std::string &transname, + const mesh &source_mesh, const std::string &source_displacements, + const mesh_region &source_region, const mesh &target_mesh, + const std::string &target_displacements, const mesh_region &target_region) + { + auto &mf_source = md.mesh_fem_of_variable(source_displacements); + auto mf_target = md.mesh_fem_of_variable(target_displacements); + auto p_transformation + = std::make_shared<interpolate_transformation_on_deformed_domains>(source_region.id(), + mf_source, + source_displacements, + target_region.id(), + mf_target, + target_displacements); + md.add_interpolate_transformation(transname, p_transformation); + } + +} /* end of namespace getfem. */ Propchange: trunk/getfem/src/getfem_interpolation_on_deformed_domains.cpp ------------------------------------------------------------------------------ svn:keywords = Id URL Header _______________________________________________ Getfem-commits mailing list Getfem-commits@gna.org https://mail.gna.org/listinfo/getfem-commits