diff --git a/src/geometry/mesh_sphere_intersector.hh b/src/geometry/mesh_sphere_intersector.hh index d9df19956..27663560e 100644 --- a/src/geometry/mesh_sphere_intersector.hh +++ b/src/geometry/mesh_sphere_intersector.hh @@ -1,101 +1,115 @@ /** * @file mesh_segment_intersector.hh * * @author Clement Roux-Langlois * * @date creation: Wed Jun 10 2015 * * @brief Computation of mesh intersection with sphere(s) * * @section LICENSE * * Copyright (©) 2010-2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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. * * Akantu 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_SPHERE_INTERSECTOR_HH__ #define __AKANTU_MESH_SPHERE_INTERSECTOR_HH__ #include "aka_common.hh" #include "mesh_geom_intersector.hh" #include "mesh_geom_common.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ +/* -------------------------------------------------------------------------- */ +/* class for new igfem elements mesh events */ +/* -------------------------------------------------------------------------- */ +#if defined(AKANTU_IGFEM) +class NewIGFEMElementsEvent : public NewElementsEvent { +public: + AKANTU_GET_MACRO_NOT_CONST(OldElementsList, old_elements, Array &); + AKANTU_GET_MACRO(OldElementsList, old_elements, const Array &); +protected: + Array old_elements; +}; +#endif + + /// Here, we know what kernel we have to use typedef Spherical SK; template class MeshSphereIntersector : public MeshGeomIntersector, SK::Sphere_3, SK> { /// Parent class type typedef MeshGeomIntersector, SK::Sphere_3, SK> parent_type; /// Result of intersection function type typedef typename IntersectionTypeHelper, K>, K::Segment_3>::intersection_type result_type; /// Pair of intersection points and element id typedef std::pair pair_type; public: /// Construct from mesh explicit MeshSphereIntersector(const Mesh & mesh); /// Destructor virtual ~MeshSphereIntersector(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the new_node_per_elem array AKANTU_GET_MACRO(NewNodePerElem, new_node_per_elem, const Array) public: /// Construct the primitive tree object virtual void constructData(); /** * @brief Computes the intersection of the mesh with a sphere * * @param query (sphere) to compute the intersections with the mesh */ virtual void computeIntersectionQuery(const SK::Sphere_3 & query); /// Build the IGFEM mesh virtual void buildResultFromQueryList(const std::list & query); /// Remove the additionnal nodes void removeAdditionnalNodes(); protected: /// new node per element Array new_node_per_elem; /// number of fem nodes in the initial mesh const UInt nb_nodes_fem; }; __END_AKANTU__ #include "mesh_sphere_intersector_tmpl.hh" #endif // __AKANTU_MESH_SPHERE_INTERSECTOR_HH__ diff --git a/src/geometry/mesh_sphere_intersector_tmpl.hh b/src/geometry/mesh_sphere_intersector_tmpl.hh index 85a7036ca..fcf0fb241 100644 --- a/src/geometry/mesh_sphere_intersector_tmpl.hh +++ b/src/geometry/mesh_sphere_intersector_tmpl.hh @@ -1,336 +1,362 @@ /** * @file mesh_sphere_intersector_tmpl.hh * * @author Clément Roux-Langlois * * @date creation: Wed june 10 2015 * @date last modification: Wed June 17 2015 * * @brief Computation of mesh intersection with spheres * * @section LICENSE * * Copyright (©) 2010-2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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. * * Akantu 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_SPHERE_INTERSECTOR_TMPL_HH__ #define __AKANTU_MESH_SPHERE_INTERSECTOR_TMPL_HH__ #include "aka_common.hh" #include "mesh_geom_common.hh" #include "tree_type_helper.hh" +#include "mesh_sphere_intersector.hh" __BEGIN_AKANTU__ template MeshSphereIntersector::MeshSphereIntersector(const Mesh & mesh): parent_type(mesh), new_node_per_elem(0, 1 + 4*(dim-1)), nb_nodes_fem(mesh.getNodes().getSize()) { this->constructData(); for(Mesh::type_iterator it = mesh.firstType(dim); it != mesh.lastType(dim); ++it){ #if defined(AKANTU_IGFEM) if(*it == _triangle_3){ // Addition of the igfem types in the mesh const_cast(this->mesh).addConnectivityType(_igfem_triangle_4, _not_ghost); const_cast(this->mesh).addConnectivityType(_igfem_triangle_4, _ghost); const_cast(this->mesh).addConnectivityType(_igfem_triangle_5, _not_ghost); const_cast(this->mesh).addConnectivityType(_igfem_triangle_5, _ghost); } else if( (*it != _triangle_3) && (*it != _igfem_triangle_4) && (*it != _igfem_triangle_5) ) AKANTU_DEBUG_ERROR("Not ready for mesh type " << *it); #else if( (*it != _triangle_3) ) AKANTU_DEBUG_ERROR("Not ready for mesh type " << *it); #endif } } template MeshSphereIntersector::~MeshSphereIntersector() {} template void MeshSphereIntersector::constructData() { this->new_node_per_elem.resize(this->mesh.getConnectivity(type).getSize()); this->new_node_per_elem.clear(); MeshGeomIntersector, SK::Sphere_3, SK>::constructData(); } template void MeshSphereIntersector::computeIntersectionQuery(const SK::Sphere_3 & query) { AKANTU_DEBUG_IN(); Array & nodes = const_cast &>(this->mesh.getNodes()); UInt nb_node = nodes.getSize() ; Real tol = 1e-10; typedef boost::variant sk_inter_res; TreeTypeHelper, Spherical>::const_iterator it = this->factory.getPrimitiveList().begin(), end= this->factory.getPrimitiveList().end(); NewNodesEvent new_nodes; for(; it!=end; ++it){ std::list s_results; CGAL::intersection(*it, query, std::back_inserter(s_results)); if(s_results.size()==1){ // just one point if (pair_type * pair = boost::get(&s_results.front())){ if(pair->second==1){ // not a point tangent to the sphere // Addition of the new node Vector new_node(dim, 0.0); SK::Circular_arc_point_3 arc_point = pair->first; new_node(0) = to_double(arc_point.x()); new_node(1) = to_double(arc_point.y()); bool is_on_mesh = false, is_new = true; UInt n = nb_nodes_fem-1; for(; n != nb_node; ++n){ if( ( pow((new_node(0)-nodes(n,0)),2) + pow((new_node(1)-nodes(n,1)),2) ) < pow(tol,2) ){ is_new = false; break; } } Real x1 = to_double(it->source().x()), y1 = to_double(it->source().y()), x2 = to_double(it->target().x()), y2 = to_double(it->target().y()), l = pow( pow((new_node(0)-x1),2) + pow((new_node(1)-y1),2), 1/2); if( ( ( pow((new_node(0)-x1),2) + pow((new_node(1)-y1),2) ) < pow(l*tol,2) ) || ( ( pow((new_node(0)-x2),2) + pow((new_node(1)-y2),2) ) < pow(l*tol,2) ) ){ is_on_mesh = true; is_new = false; } if(is_new){ nodes.push_back(new_node); new_nodes.getList().push_back(nb_node); nb_node += 1 ; } if(!is_on_mesh){ new_node_per_elem(it->id(), 0) += 1; new_node_per_elem(it->id(), ( 2*new_node_per_elem(it->id(),0)) - 1 ) = n ; new_node_per_elem(it->id(), 2*new_node_per_elem(it->id(),0) ) = it->segId() ; } } } } } - const_cast(this->mesh).sendEvent(new_nodes); + if(new_nodes.getList().getSize()) + const_cast(this->mesh).sendEvent(new_nodes); AKANTU_DEBUG_OUT(); } template void MeshSphereIntersector::removeAdditionnalNodes() { AKANTU_DEBUG_IN(); RemovedNodesEvent remove_nodes(this->mesh); Array & nodes_removed = remove_nodes.getList(); Array & new_numbering = remove_nodes.getNewNumbering(); UInt nnod = 0; for(; nnod != nb_nodes_fem ; ++nnod){ new_numbering(nnod) = nnod ; } for(nnod = nb_nodes_fem ; nnod != this->mesh.getNodes().getSize(); ++nnod){ new_numbering(nnod) = UInt(-1) ; nodes_removed.push_back(nnod); } const_cast(this->mesh).sendEvent(remove_nodes); AKANTU_DEBUG_OUT(); } #if defined(AKANTU_IGFEM) template void MeshSphereIntersector::buildResultFromQueryList(const std::list & query_list) { AKANTU_DEBUG_IN(); this->computeIntersectionQueryList(query_list); Array connec_type_tmpl = this->mesh.getConnectivity(type); Array & connec_igfem_tri4 = const_cast &>(this->mesh.getConnectivity(_igfem_triangle_4)), & connec_igfem_tri5 = const_cast &>(this->mesh.getConnectivity(_igfem_triangle_5)), & connec_tri3 = const_cast &>(this->mesh.getConnectivity(_triangle_3)); Element element_tri3(_triangle_3, 0, _not_ghost, _ek_regular), element_tri4(_igfem_triangle_4, 0, _not_ghost, _ek_igfem), element_tri5(_igfem_triangle_5, 0, _not_ghost, _ek_igfem); - NewElementsEvent new_elements; + NewIGFEMElementsEvent new_elements; UInt n_new_el = 0; - Array & new_elements_list = new_elements.getList() ; - new_elements.getList().extendComponentsInterlaced(2, 1); + Array & new_elements_list = new_elements.getList(); + Array & old_elements_list = new_elements.getOldElementsList(); + RemovedElementsEvent remove_elem(this->mesh); remove_elem.getNewNumbering().alloc(connec_type_tmpl.getSize(), 1, type, _not_ghost); Array & new_numbering = remove_elem.getNewNumbering(type, _not_ghost); UInt new_nb_type_tmpl = 0; // Meter of the number of element (type) will we loop on elements for(UInt nel=0; nel != new_node_per_elem.getSize(); ++nel){ if( (type != _triangle_3) && (new_node_per_elem(nel,0)==0) ){ - Element element_type_tmpl(type, 0, _not_ghost); + /// this element was previously an IGFEM element and is replaced by a regular element + /// store the old element that will be removed + Element element_type_tmpl(type, 0, _not_ghost, Mesh::getKind(type)); new_elements_list.resize(n_new_el+1); + old_elements_list.resize(n_new_el+1); Vector ctri3(3); ctri3(0) = connec_type_tmpl(nel,0); ctri3(1) = connec_type_tmpl(nel,1); ctri3(2) = connec_type_tmpl(nel,2); + /// add the new element in the mesh UInt nb_tri3 = connec_tri3.getSize(); connec_tri3.push_back(ctri3); element_tri3.element = nb_tri3; - new_elements_list(n_new_el,0) = element_tri3; + new_elements_list(n_new_el) = element_tri3; + /// the number of the old element in the mesh element_type_tmpl.element = nel; - new_elements_list(n_new_el,1) = element_type_tmpl; + old_elements_list(n_new_el) = element_type_tmpl; new_numbering(nel) = UInt(-1); ++n_new_el; + + remove_elem.getList().push_back(element_type_tmpl); + } else if(new_node_per_elem(nel,0)!=0){ Element element_type_tmpl(type, 0, _not_ghost); + element_type_tmpl.kind = Mesh::getKind(type); new_elements_list.resize(n_new_el+1); + old_elements_list.resize(n_new_el+1); switch(new_node_per_elem(nel,0)){ case 1 :{ Vector ctri4(4); switch(new_node_per_elem(nel,2)){ case 1 : ctri4(0) = connec_type_tmpl(nel,2); ctri4(1) = connec_type_tmpl(nel,0); ctri4(2) = connec_type_tmpl(nel,1); break; case 2 : ctri4(0) = connec_type_tmpl(nel,0); ctri4(1) = connec_type_tmpl(nel,1); ctri4(2) = connec_type_tmpl(nel,2); break; case 3 : ctri4(0) = connec_type_tmpl(nel,1); ctri4(1) = connec_type_tmpl(nel,2); ctri4(2) = connec_type_tmpl(nel,0); break; default : AKANTU_DEBUG_ERROR("A triangle have only 3 segments not "<< new_node_per_elem(nel,2)); break; } ctri4(3) = new_node_per_elem(nel,1); UInt nb_tri4 = connec_igfem_tri4.getSize(); connec_igfem_tri4.push_back(ctri4); element_tri4.element = nb_tri4; //new_elements.getList().push_back(element_tri4); - new_elements_list(n_new_el,0) = element_tri4; + new_elements_list(n_new_el) = element_tri4; if(type == _igfem_triangle_4) new_numbering.push_back(connec_igfem_tri4.getSize()-2); break; } case 2 :{ Vector ctri5(5); if( (new_node_per_elem(nel,2)==1) && (new_node_per_elem(nel,4)==2) ){ ctri5(0) = connec_type_tmpl(nel,1); ctri5(1) = connec_type_tmpl(nel,2); ctri5(2) = connec_type_tmpl(nel,0); ctri5(3) = new_node_per_elem(nel,3); ctri5(4) = new_node_per_elem(nel,1); } else if((new_node_per_elem(nel,2)==1) && (new_node_per_elem(nel,4)==3)){ ctri5(0) = connec_type_tmpl(nel,0); ctri5(1) = connec_type_tmpl(nel,1); ctri5(2) = connec_type_tmpl(nel,2); ctri5(3) = new_node_per_elem(nel,1); ctri5(4) = new_node_per_elem(nel,3); } else if((new_node_per_elem(nel,2)==2) && (new_node_per_elem(nel,4)==3)){ ctri5(0) = connec_type_tmpl(nel,2); ctri5(1) = connec_type_tmpl(nel,0); ctri5(2) = connec_type_tmpl(nel,1); ctri5(3) = new_node_per_elem(nel,3); ctri5(4) = new_node_per_elem(nel,1); } else if((new_node_per_elem(nel,2)==2) && (new_node_per_elem(nel,4)==1)){ ctri5(0) = connec_type_tmpl(nel,1); ctri5(1) = connec_type_tmpl(nel,2); ctri5(2) = connec_type_tmpl(nel,0); ctri5(3) = new_node_per_elem(nel,1); ctri5(4) = new_node_per_elem(nel,3); } else if((new_node_per_elem(nel,2)==3) && (new_node_per_elem(nel,4)==1)){ ctri5(0) = connec_type_tmpl(nel,0); ctri5(1) = connec_type_tmpl(nel,1); ctri5(2) = connec_type_tmpl(nel,2); ctri5(3) = new_node_per_elem(nel,3); ctri5(4) = new_node_per_elem(nel,1); } else if((new_node_per_elem(nel,2)==3) && (new_node_per_elem(nel,4)==2)){ ctri5(0) = connec_type_tmpl(nel,2); ctri5(1) = connec_type_tmpl(nel,0); ctri5(2) = connec_type_tmpl(nel,1); ctri5(3) = new_node_per_elem(nel,1); ctri5(4) = new_node_per_elem(nel,3); } else{ AKANTU_DEBUG_ERROR("A triangle have only 3 segments not "<< new_node_per_elem(nel,2) << "and" << new_node_per_elem(nel,4)); } - connec_igfem_tri5.push_back(ctri5); UInt nb_tri5 = connec_igfem_tri5.getSize(); + connec_igfem_tri5.push_back(ctri5); element_tri5.element = nb_tri5; //new_elements.getList().push_back(element_tri5); - new_elements_list(n_new_el,0) = element_tri5; + new_elements_list(n_new_el) = element_tri5; if(type == _igfem_triangle_5){ new_numbering.push_back(new_nb_type_tmpl); ++new_nb_type_tmpl; } break; } default: AKANTU_DEBUG_ERROR("Igfem cannot add "<< new_node_per_elem(nel,0) << " nodes to triangles"); break; } element_type_tmpl.element = nel; - new_elements_list(n_new_el,1) = element_type_tmpl; + old_elements_list(n_new_el) = element_type_tmpl; remove_elem.getList().push_back(element_type_tmpl); new_numbering(nel) = UInt(-1); ++n_new_el; } else{ new_numbering(nel) = new_nb_type_tmpl; ++new_nb_type_tmpl; } } + // for(UInt nel=new_node_per_elem.getSize(); nel < this->mesh.getNbElement(type); ++nel) { + // new_numbering(nel) = new_nb_type_tmpl; + // ++new_nb_type_tmpl; + // } + + UInt el_index = 0; + for (UInt e = 0; e < this->mesh.getNbElement(type); ++e) { + if (new_numbering(e) != UInt(-1)) { + new_numbering(e) = el_index; + ++el_index; + } + } + if(n_new_el > 0){ const_cast(this->mesh).sendEvent(new_elements); const_cast(this->mesh).sendEvent(remove_elem); } AKANTU_DEBUG_OUT(); } #else template void MeshSphereIntersector::buildResultFromQueryList(const std::list & query_list) { AKANTU_DEBUG_TO_IMPLEMENT(); } #endif __END_AKANTU__ #endif // __AKANTU_MESH_SPHERE_INTERSECTOR_TMPL_HH__ diff --git a/src/model/solid_mechanics/material.cc b/src/model/solid_mechanics/material.cc index 1f872e4bc..7431c8159 100644 --- a/src/model/solid_mechanics/material.cc +++ b/src/model/solid_mechanics/material.cc @@ -1,1803 +1,1872 @@ /** * @file material.cc * * @author Aurelia Isabel Cuba Ramos * @author Marco Vocialta * @author Nicolas Richart * @author Daniel Pino Muñoz * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Sep 16 2014 * * @brief Implementation of the common part of the material class * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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. * * Akantu 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "solid_mechanics_model.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, const ID & id) : Memory(id, model.getMemoryID()), Parsable(_st_material, id), is_init(false), fem(&(model.getFEEngine())), finite_deformation(false), name(""), model(&model), spatial_dimension(this->model->getSpatialDimension()), element_filter("element_filter", id, this->memory_id), stress("stress", *this), eigenstrain("eigenstrain", *this), gradu("grad_u", *this), green_strain("green_strain",*this), piola_kirchhoff_2("piola_kirchhoff_2", *this), // potential_energy_vector(false), potential_energy("potential_energy", *this), is_non_local(false), use_previous_stress(false), use_previous_gradu(false), interpolation_inverse_coordinates("interpolation inverse coordinates", *this), interpolation_points_matrices("interpolation points matrices", *this) { AKANTU_DEBUG_IN(); /// for each connectivity types allocate the element filer array of the material model.getMesh().initElementTypeMapArray(element_filter, 1, spatial_dimension, false, _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id) : Memory(id, model.getMemoryID()), Parsable(_st_material, id), is_init(false), fem(&(model.getFEEngine())), finite_deformation(false), name(""), model(&model), spatial_dimension(dim), element_filter("element_filter", id, this->memory_id), stress("stress", *this, dim, fe_engine, this->element_filter), eigenstrain("eignestrain", *this, dim, fe_engine, this->element_filter), gradu("gradu", *this, dim, fe_engine, this->element_filter), green_strain("green_strain", *this, dim, fe_engine, this->element_filter), piola_kirchhoff_2("poila_kirchhoff_2", *this, dim, fe_engine, this->element_filter), potential_energy("potential_energy", *this, dim, fe_engine, this->element_filter), is_non_local(false), use_previous_stress(false), use_previous_gradu(false), interpolation_inverse_coordinates("interpolation inverse_coordinates", *this, dim, fe_engine, this->element_filter), interpolation_points_matrices("interpolation points matrices", *this, dim, fe_engine, this->element_filter) { AKANTU_DEBUG_IN(); mesh.initElementTypeMapArray(element_filter, 1, spatial_dimension, false, _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::~Material() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::initialize() { registerParam("rho" , rho , 0. , _pat_parsable | _pat_modifiable, "Density"); registerParam("name" , name , std::string(), _pat_parsable | _pat_readable); registerParam("finite_deformation" , finite_deformation , false , _pat_parsable | _pat_readable, "Is finite deformation"); registerParam("inelastic_deformation", inelastic_deformation, false , _pat_internal, "Is inelastic deformation"); /// allocate gradu stress for local elements eigenstrain.initialize(spatial_dimension * spatial_dimension); gradu.initialize(spatial_dimension * spatial_dimension); stress.initialize(spatial_dimension * spatial_dimension); this->model->registerEventHandler(*this); } /* -------------------------------------------------------------------------- */ void Material::initMaterial() { AKANTU_DEBUG_IN(); if(finite_deformation) { this->piola_kirchhoff_2.initialize(spatial_dimension * spatial_dimension); if(use_previous_stress) this->piola_kirchhoff_2.initializeHistory(); this->green_strain.initialize(spatial_dimension * spatial_dimension); } if(use_previous_stress) this->stress.initializeHistory(); if(use_previous_gradu) this->gradu.initializeHistory(); for (std::map *>::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->resize(); for (std::map *>::iterator it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->resize(); for (std::map *>::iterator it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->resize(); is_init = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::savePreviousState() { AKANTU_DEBUG_IN(); for (std::map *>::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) { if(it->second->hasHistory()) it->second->saveCurrentValues(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the residual by assembling @f$\int_{e} \sigma_e \frac{\partial * \varphi}{\partial X} dX @f$ * * @param[in] displacements nodes displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::updateResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); computeAllStresses(ghost_type); assembleResidual(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::assembleResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); if(!finite_deformation){ Array & residual = const_cast &>(model->getResidual()); Mesh & mesh = fem->getMesh(); Mesh::type_iterator it = element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element) { const Array & shapes_derivatives = fem->getShapesDerivatives(*it, ghost_type); UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt nb_quadrature_points = fem->getNbQuadraturePoints(*it, ghost_type); /// compute @f$\sigma \frac{\partial \varphi}{\partial X}@f$ by @f$\mathbf{B}^t \mathbf{\sigma}_q@f$ Array * sigma_dphi_dx = new Array(nb_element*nb_quadrature_points, size_of_shapes_derivatives, "sigma_x_dphi_/_dX"); Array * shapesd_filtered = new Array(0, size_of_shapes_derivatives, "filtered shapesd"); FEEngine::filterElementalData(mesh, shapes_derivatives, *shapesd_filtered, *it, ghost_type, elem_filter); Array & stress_vect = this->stress(*it, ghost_type); Array::matrix_iterator sigma = stress_vect.begin(spatial_dimension, spatial_dimension); Array::matrix_iterator B = shapesd_filtered->begin(spatial_dimension, nb_nodes_per_element); Array::matrix_iterator Bt_sigma_it = sigma_dphi_dx->begin(spatial_dimension, nb_nodes_per_element); for (UInt q = 0; q < nb_element*nb_quadrature_points; ++q, ++sigma, ++B, ++Bt_sigma_it) Bt_sigma_it->mul(*sigma, *B); delete shapesd_filtered; /** * compute @f$\int \sigma * \frac{\partial \varphi}{\partial X}dX@f$ by @f$ \sum_q \mathbf{B}^t * \mathbf{\sigma}_q \overline w_q J_q@f$ */ Array * int_sigma_dphi_dx = new Array(nb_element, nb_nodes_per_element * spatial_dimension, "int_sigma_x_dphi_/_dX"); fem->integrate(*sigma_dphi_dx, *int_sigma_dphi_dx, size_of_shapes_derivatives, *it, ghost_type, elem_filter); delete sigma_dphi_dx; /// assemble fem->assembleArray(*int_sigma_dphi_dx, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), residual.getNbComponent(), *it, ghost_type, elem_filter, -1); delete int_sigma_dphi_dx; } } } else{ switch (spatial_dimension){ case 1: this->assembleResidual<1>(ghost_type); break; case 2: this->assembleResidual<2>(ghost_type); break; case 3: this->assembleResidual<3>(ghost_type); break; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stress from the gradu * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::computeAllStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Mesh::type_iterator it = fem->getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = fem->getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); if (elem_filter.getSize() == 0) continue; Array & gradu_vect = gradu(*it, ghost_type); /// compute @f$\nabla u@f$ fem->gradientOnQuadraturePoints(model->getDisplacement(), gradu_vect, spatial_dimension, *it, ghost_type, elem_filter); gradu_vect -= eigenstrain(*it, ghost_type); /// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$ computeStress(*it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computeAllCauchyStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(finite_deformation,"The Cauchy stress can only be computed if you are working in finite deformation."); //resizeInternalArray(stress); Mesh::type_iterator it = fem->getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = fem->getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) switch (spatial_dimension){ case 1: this->computeCauchyStress<1>(*it, ghost_type); break; case 2: this->computeCauchyStress<2>(*it, ghost_type); break; case 3: this->computeCauchyStress<3>(*it, ghost_type); break; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::computeCauchyStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array::matrix_iterator gradu_it = this->gradu(el_type, ghost_type).begin(dim, dim); Array::matrix_iterator gradu_end = this->gradu(el_type, ghost_type).end(dim, dim); Array::matrix_iterator piola_it = this->piola_kirchhoff_2(el_type, ghost_type).begin(dim, dim); Array::matrix_iterator stress_it = this->stress(el_type, ghost_type).begin(dim, dim); Matrix F_tensor(dim, dim); for (; gradu_it != gradu_end; ++gradu_it, ++piola_it, ++stress_it) { Matrix & grad_u = *gradu_it; Matrix & piola = *piola_it; Matrix & sigma = *stress_it; gradUToF (grad_u, F_tensor); this->computeCauchyStressOnQuad(F_tensor, piola, sigma); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::setToSteadyState(GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & displacement = model->getDisplacement(); //resizeInternalArray(gradu); UInt spatial_dimension = model->getSpatialDimension(); Mesh::type_iterator it = fem->getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = fem->getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); Array & gradu_vect = gradu(*it, ghost_type); /// compute @f$\nabla u@f$ fem->gradientOnQuadraturePoints(displacement, gradu_vect, spatial_dimension, *it, ghost_type, elem_filter); setToSteadyState(*it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stiffness matrix by assembling @f$\int_{\omega} B^t \times D * \times B d\omega @f$ * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Mesh::type_iterator it = element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { if(finite_deformation){ switch (spatial_dimension) { case 1: { assembleStiffnessMatrixNL < 1 > (*it, ghost_type); assembleStiffnessMatrixL2 < 1 > (*it, ghost_type); break; } case 2: { assembleStiffnessMatrixNL < 2 > (*it, ghost_type); assembleStiffnessMatrixL2 < 2 > (*it, ghost_type); break; } case 3: { assembleStiffnessMatrixNL < 3 > (*it, ghost_type); assembleStiffnessMatrixL2 < 3 > (*it, ghost_type); break; } } } else { switch(spatial_dimension) { case 1: { assembleStiffnessMatrix<1>(*it, ghost_type); break; } case 2: { assembleStiffnessMatrix<2>(*it, ghost_type); break; } case 3: { assembleStiffnessMatrix<3>(*it, ghost_type); break; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.getSize()) { SparseMatrix & K = const_cast(model->getStiffnessMatrix()); const Array & shapes_derivatives = fem->getShapesDerivatives(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem->getNbQuadraturePoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); fem->gradientOnQuadraturePoints(model->getDisplacement(), gradu_vect, dim, type, ghost_type, elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); Array * tangent_stiffness_matrix = new Array(nb_element*nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); tangent_stiffness_matrix->clear(); computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type); Array * shapesd_filtered = new Array(0, dim * nb_nodes_per_element, "filtered shapesd"); FEEngine::filterElementalData(fem->getMesh(), shapes_derivatives, *shapesd_filtered, type, ghost_type, elem_filter); /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; Array * bt_d_b = new Array(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); Matrix B(tangent_size, dim * nb_nodes_per_element); Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); Array::matrix_iterator shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::matrix_iterator Bt_D_B_it = bt_d_b->begin(dim*nb_nodes_per_element, dim*nb_nodes_per_element); Array::matrix_iterator D_it = tangent_stiffness_matrix->begin(tangent_size, tangent_size); Array::matrix_iterator D_end = tangent_stiffness_matrix->end (tangent_size, tangent_size); for(; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it) { Matrix & D = *D_it; Matrix & Bt_D_B = *Bt_D_B_it; VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bt_D.mul(B, D); Bt_D_B.mul(Bt_D, B); } delete tangent_stiffness_matrix; delete shapesd_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array(nb_element, bt_d_b_size * bt_d_b_size, "K_e"); fem->integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, elem_filter); delete bt_d_b; fem->assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, elem_filter); delete K_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixNL(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); SparseMatrix & K = const_cast (model->getStiffnessMatrix()); const Array & shapes_derivatives = fem->getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); //Array & gradu_vect = delta_gradu(type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem->getNbQuadraturePoints(type, ghost_type); //gradu_vect.resize(nb_quadrature_points * nb_element); // fem->gradientOnQuadraturePoints(model->getIncrement(), gradu_vect, // dim, type, ghost_type, &elem_filter); Array * shapes_derivatives_filtered = new Array (nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); Array::const_matrix_iterator shapes_derivatives_it = shapes_derivatives.begin(spatial_dimension, nb_nodes_per_element); Array::matrix_iterator shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element); UInt * elem_filter_val = elem_filter.storage(); for (UInt e = 0; e < nb_element; ++e, ++elem_filter_val) for (UInt q = 0; q < nb_quadrature_points; ++q, ++shapes_derivatives_filtered_it) *shapes_derivatives_filtered_it = shapes_derivatives_it[*elem_filter_val * nb_quadrature_points + q]; /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_s_b_size = dim * nb_nodes_per_element; Array * bt_s_b = new Array (nb_element * nb_quadrature_points, bt_s_b_size * bt_s_b_size, "B^t*D*B"); UInt piola_matrix_size = getCauchyStressMatrixSize(dim); Matrix B(piola_matrix_size, bt_s_b_size); Matrix Bt_S(bt_s_b_size, piola_matrix_size); Matrix S(piola_matrix_size, piola_matrix_size); shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element); Array::matrix_iterator Bt_S_B_it = bt_s_b->begin(bt_s_b_size, bt_s_b_size); Array::matrix_iterator Bt_S_B_end = bt_s_b->end(bt_s_b_size, bt_s_b_size); Array::matrix_iterator piola_it = piola_kirchhoff_2(type, ghost_type).begin(dim, dim); for (; Bt_S_B_it != Bt_S_B_end; ++Bt_S_B_it, ++shapes_derivatives_filtered_it, ++piola_it) { Matrix & Bt_S_B = *Bt_S_B_it; Matrix & Piola_kirchhoff_matrix = *piola_it; setCauchyStressMatrix< dim >(Piola_kirchhoff_matrix, S); VoigtHelper::transferBMatrixToBNL(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bt_S.mul < true, false > (B, S); Bt_S_B.mul < false, false > (Bt_S, B); } delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array (nb_element, bt_s_b_size * bt_s_b_size, "K_e"); fem->integrate(*bt_s_b, *K_e, bt_s_b_size * bt_s_b_size, type, ghost_type, elem_filter); delete bt_s_b; fem->assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixL2(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); SparseMatrix & K = const_cast (model->getStiffnessMatrix()); const Array & shapes_derivatives = fem->getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem->getNbQuadraturePoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); fem->gradientOnQuadraturePoints(model->getDisplacement(), gradu_vect, dim, type, ghost_type, elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); Array * tangent_stiffness_matrix = new Array (nb_element*nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); tangent_stiffness_matrix->clear(); computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type); Array * shapes_derivatives_filtered = new Array (nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); Array::const_matrix_iterator shapes_derivatives_it = shapes_derivatives.begin(spatial_dimension, nb_nodes_per_element); Array::matrix_iterator shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element); UInt * elem_filter_val = elem_filter.storage(); for (UInt e = 0; e < nb_element; ++e, ++elem_filter_val) for (UInt q = 0; q < nb_quadrature_points; ++q, ++shapes_derivatives_filtered_it) *shapes_derivatives_filtered_it = shapes_derivatives_it[*elem_filter_val * nb_quadrature_points + q]; /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; Array * bt_d_b = new Array (nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); Matrix B(tangent_size, dim * nb_nodes_per_element); Matrix B2(tangent_size, dim * nb_nodes_per_element); Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element); Array::matrix_iterator Bt_D_B_it = bt_d_b->begin(dim*nb_nodes_per_element, dim * nb_nodes_per_element); Array::matrix_iterator grad_u_it = gradu_vect.begin(dim, dim); Array::matrix_iterator D_it = tangent_stiffness_matrix->begin(tangent_size, tangent_size); Array::matrix_iterator D_end = tangent_stiffness_matrix->end(tangent_size, tangent_size); for (; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it, ++grad_u_it) { Matrix & grad_u = *grad_u_it; Matrix & D = *D_it; Matrix & Bt_D_B = *Bt_D_B_it; //transferBMatrixToBL1 (*shapes_derivatives_filtered_it, B, nb_nodes_per_element); VoigtHelper::transferBMatrixToSymVoigtBMatrix(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); VoigtHelper::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2, nb_nodes_per_element); B += B2; Bt_D.mul < true, false > (B, D); Bt_D_B.mul < false, false > (Bt_D, B); } delete tangent_stiffness_matrix; delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array (nb_element, bt_d_b_size * bt_d_b_size, "K_e"); fem->integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, elem_filter); delete bt_d_b; fem->assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleResidual(GhostType ghost_type){ AKANTU_DEBUG_IN(); Array & residual = const_cast &> (model->getResidual()); Mesh & mesh = fem->getMesh(); Mesh::type_iterator it = element_filter.firstType(dim, ghost_type); Mesh::type_iterator last_type = element_filter.lastType(dim, ghost_type); for (; it != last_type; ++it) { const Array & shapes_derivatives = fem->getShapesDerivatives(*it, ghost_type); Array & elem_filter = element_filter(*it, ghost_type); if (elem_filter.getSize() == 0) continue; UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt nb_quadrature_points = fem->getNbQuadraturePoints(*it, ghost_type); Array * shapesd_filtered = new Array(0, size_of_shapes_derivatives, "filtered shapesd"); FEEngine::filterElementalData(mesh, shapes_derivatives, *shapesd_filtered, *it, ghost_type, elem_filter); Array::matrix_iterator shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); //Set stress vectors UInt stress_size = getTangentStiffnessVoigtSize(dim); //Set matrices B and BNL* UInt bt_s_size = dim * nb_nodes_per_element; Array * bt_s = new Array (nb_element * nb_quadrature_points, bt_s_size, "B^t*S"); Array::matrix_iterator grad_u_it = this->gradu(*it, ghost_type).begin(dim, dim); Array::matrix_iterator grad_u_end = this->gradu(*it, ghost_type).end(dim, dim); Array::matrix_iterator stress_it = this->piola_kirchhoff_2(*it, ghost_type).begin(dim, dim); shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::matrix_iterator bt_s_it = bt_s->begin(bt_s_size, 1); Matrix S_vect(stress_size, 1); Matrix B_tensor(stress_size, bt_s_size); Matrix B2_tensor(stress_size, bt_s_size); for (; grad_u_it != grad_u_end; ++grad_u_it, ++stress_it, ++shapes_derivatives_filtered_it, ++bt_s_it) { Matrix & grad_u = *grad_u_it; Matrix & r_it = *bt_s_it; Matrix & S_it = *stress_it; SetCauchyStressArray (S_it, S_vect); VoigtHelper::transferBMatrixToSymVoigtBMatrix(*shapes_derivatives_filtered_it, B_tensor, nb_nodes_per_element); VoigtHelper::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2_tensor, nb_nodes_per_element); B_tensor += B2_tensor; r_it.mul < true, false > (B_tensor, S_vect); } delete shapesd_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * r_e = new Array (nb_element, bt_s_size, "r_e"); fem->integrate(*bt_s, *r_e, bt_s_size, *it, ghost_type, elem_filter); delete bt_s; fem->assembleArray(*r_e, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), residual.getNbComponent(), *it, ghost_type, elem_filter, -1); delete r_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computeAllStressesFromTangentModuli(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Mesh::type_iterator it = element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { switch(spatial_dimension) { case 1: { computeAllStressesFromTangentModuli<1>(*it, ghost_type); break; } case 2: { computeAllStressesFromTangentModuli<2>(*it, ghost_type); break; } case 3: { computeAllStressesFromTangentModuli<3>(*it, ghost_type); break; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::computeAllStressesFromTangentModuli(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & shapes_derivatives = fem->getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem->getNbQuadraturePoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); Array & disp = model->getDisplacement(); fem->gradientOnQuadraturePoints(disp, gradu_vect, dim, type, ghost_type, elem_filter); UInt tangent_moduli_size = getTangentStiffnessVoigtSize(dim); Array * tangent_moduli_tensors = new Array(nb_element*nb_quadrature_points, tangent_moduli_size * tangent_moduli_size, "tangent_moduli_tensors"); tangent_moduli_tensors->clear(); computeTangentModuli(type, *tangent_moduli_tensors, ghost_type); Array * shapesd_filtered = new Array(0, dim* nb_nodes_per_element, "filtered shapesd"); FEEngine::filterElementalData(fem->getMesh(), shapes_derivatives, *shapesd_filtered, type, ghost_type, elem_filter); Array filtered_u(nb_element, nb_nodes_per_element * spatial_dimension); FEEngine::extractNodalToElementField(fem->getMesh(), disp, filtered_u, type, ghost_type, elem_filter); /// compute @f$\mathbf{D} \mathbf{B} \mathbf{u}@f$ Array::matrix_iterator shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::matrix_iterator D_it = tangent_moduli_tensors->begin(tangent_moduli_size, tangent_moduli_size); Array::matrix_iterator sigma_it = stress(type, ghost_type).begin(spatial_dimension, spatial_dimension); Array::vector_iterator u_it = filtered_u.begin(spatial_dimension * nb_nodes_per_element); Matrix B(tangent_moduli_size, spatial_dimension * nb_nodes_per_element); Vector Bu(tangent_moduli_size); Vector DBu(tangent_moduli_size); for (UInt e = 0; e < nb_element; ++e, ++u_it) { for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it, ++shapes_derivatives_filtered_it, ++sigma_it) { Vector & u = *u_it; Matrix & sigma = *sigma_it; Matrix & D = *D_it; VoigtHelper::transferBMatrixToSymVoigtBMatrix(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bu.mul(B, u); DBu.mul(D, Bu); // Voigt notation to full symmetric tensor for (UInt i = 0; i < dim; ++i) sigma(i, i) = DBu(i); if(dim == 2) { sigma(0,1) = sigma(1,0) = DBu(2); } else if(dim == 3) { sigma(1,2) = sigma(2,1) = DBu(3); sigma(0,2) = sigma(2,0) = DBu(4); sigma(0,1) = sigma(1,0) = DBu(5); } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergyByElements() { AKANTU_DEBUG_IN(); Mesh::type_iterator it = element_filter.firstType(spatial_dimension); Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension); for(; it != last_type; ++it) { computePotentialEnergy(*it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergy(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if(!potential_energy.exists(el_type, ghost_type)) { UInt nb_element = element_filter(el_type, ghost_type).getSize(); UInt nb_quadrature_points = fem->getNbQuadraturePoints(el_type, _not_ghost); potential_energy.alloc(nb_element * nb_quadrature_points, 1, el_type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real epot = 0.; computePotentialEnergyByElements(); /// integrate the potential energy for each type of elements Mesh::type_iterator it = element_filter.firstType(spatial_dimension); Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension); for(; it != last_type; ++it) { epot += fem->integrate(potential_energy(*it, _not_ghost), *it, _not_ghost, element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy(ElementType & type, UInt index) { AKANTU_DEBUG_IN(); Real epot = 0.; Vector epot_on_quad_points(fem->getNbQuadraturePoints(type)); computePotentialEnergyByElement(type, index, epot_on_quad_points); epot = fem->integrate(epot_on_quad_points, type, element_filter(type)(index)); AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(std::string type) { AKANTU_DEBUG_IN(); if(type == "potential") return getPotentialEnergy(); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(std::string energy_id, ElementType type, UInt index) { AKANTU_DEBUG_IN(); if(energy_id == "potential") return getPotentialEnergy(type, index); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ void Material::computeQuadraturePointsCoordinates(ElementTypeMapArray & quadrature_points_coordinates, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); const Mesh & mesh = this->model->getMesh(); Array nodes_coordinates(mesh.getNodes(), true); nodes_coordinates += this->model->getDisplacement(); Mesh::type_iterator it = this->element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = this->element_filter.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { const Array & elem_filter = this->element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; UInt nb_tot_quad = this->fem->getNbQuadraturePoints(*it, ghost_type) * nb_element; Array & quads = quadrature_points_coordinates(*it, ghost_type); quads.resize(nb_tot_quad); this->fem->interpolateOnQuadraturePoints(nodes_coordinates, quads, spatial_dimension, *it, ghost_type, elem_filter); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::initElementalFieldInterpolation(const ElementTypeMapArray & interpolation_points_coordinates) { AKANTU_DEBUG_IN(); const Mesh & mesh = fem->getMesh(); ElementTypeMapArray quadrature_points_coordinates("quadrature_points_coordinates_for_interpolation", getID()); mesh.initElementTypeMapArray(quadrature_points_coordinates, spatial_dimension, spatial_dimension); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; computeQuadraturePointsCoordinates(quadrature_points_coordinates, ghost_type); Mesh::type_iterator it = element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last = element_filter.lastType(spatial_dimension, ghost_type); for (; it != last; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type, ghost_type); if (nb_element == 0) continue; const Array & interp_points_coord = interpolation_points_coordinates(type, ghost_type); UInt nb_interpolation_points_per_elem = interp_points_coord.getSize() / nb_element; AKANTU_DEBUG_ASSERT(interp_points_coord.getSize() % nb_element == 0, "Number of interpolation points is wrong"); #define AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD(type) \ initElementalFieldInterpolation(quadrature_points_coordinates(type, ghost_type), \ interp_points_coord, \ nb_interpolation_points_per_elem, \ ghost_type) \ AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD); #undef AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::initElementalFieldInterpolation(const Array & quad_coordinates, const Array & interpolation_points_coordinates, const UInt nb_interpolation_points_per_elem, const GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & elem_fil = element_filter(type, ghost_type); UInt nb_element = elem_fil.getSize(); UInt nb_quad_per_element = fem->getNbQuadraturePoints(type, ghost_type); if(!interpolation_inverse_coordinates.exists(type, ghost_type)) interpolation_inverse_coordinates.alloc(nb_element, nb_quad_per_element*nb_quad_per_element, type, ghost_type); else interpolation_inverse_coordinates(type, ghost_type).resize(nb_element); if(!interpolation_points_matrices.exists(type, ghost_type)) interpolation_points_matrices.alloc(nb_element, nb_interpolation_points_per_elem * nb_quad_per_element, type, ghost_type); else interpolation_points_matrices(type, ghost_type).resize(nb_element); Array & interp_inv_coord = interpolation_inverse_coordinates(type, ghost_type); Array & interp_points_mat = interpolation_points_matrices(type, ghost_type); Matrix quad_coord_matrix(nb_quad_per_element, nb_quad_per_element); Array::const_matrix_iterator quad_coords_it = quad_coordinates.begin_reinterpret(spatial_dimension, nb_quad_per_element, nb_element); Array::const_matrix_iterator points_coords_begin = interpolation_points_coordinates.begin_reinterpret(spatial_dimension, nb_interpolation_points_per_elem, interpolation_points_coordinates.getSize() / nb_interpolation_points_per_elem); Array::matrix_iterator inv_quad_coord_it = interp_inv_coord.begin(nb_quad_per_element, nb_quad_per_element); Array::matrix_iterator inv_points_mat_it = interp_points_mat.begin(nb_interpolation_points_per_elem, nb_quad_per_element); /// loop over the elements of the current material and element type for (UInt el = 0; el < nb_element; ++el, ++inv_quad_coord_it, ++inv_points_mat_it, ++quad_coords_it) { /// matrix containing the quadrature points coordinates const Matrix & quad_coords = *quad_coords_it; /// matrix to store the matrix inversion result Matrix & inv_quad_coord_matrix = *inv_quad_coord_it; /// insert the quad coordinates in a matrix compatible with the interpolation buildElementalFieldInterpolationCoodinates(quad_coords, quad_coord_matrix); /// invert the interpolation matrix inv_quad_coord_matrix.inverse(quad_coord_matrix); /// matrix containing the interpolation points coordinates const Matrix & points_coords = points_coords_begin[elem_fil(el)]; /// matrix to store the interpolation points coordinates /// compatible with these functions Matrix & inv_points_coord_matrix = *inv_points_mat_it; /// insert the quad coordinates in a matrix compatible with the interpolation buildElementalFieldInterpolationCoodinates(points_coords, inv_points_coord_matrix); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::interpolateStress(ElementTypeMapArray & result, const GhostType ghost_type) { interpolateElementalField(stress, result, ghost_type); } /* -------------------------------------------------------------------------- */ void Material::interpolateElementalField(const ElementTypeMapArray & field, ElementTypeMapArray & result, const GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh::type_iterator it = element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last = element_filter.lastType(spatial_dimension, ghost_type); for (; it != last; ++it) { ElementType type = *it; Array & elem_fil = element_filter(type, ghost_type); UInt nb_element = elem_fil.getSize(); UInt nb_quad_per_element = fem->getNbQuadraturePoints(type, ghost_type); const Array & field_vec = field(type, ghost_type); Array & result_vec = result(type, ghost_type); Matrix coefficients(nb_quad_per_element, field_vec.getNbComponent()); const Array & interp_inv_coord = interpolation_inverse_coordinates(type, ghost_type); const Array & interp_points_coord = interpolation_points_matrices(type, ghost_type); UInt nb_interpolation_points_per_elem = interp_points_coord.getNbComponent() / nb_quad_per_element; Array::const_matrix_iterator field_it = field_vec.begin_reinterpret(field_vec.getNbComponent(), nb_quad_per_element, nb_element); Array::const_matrix_iterator interpolation_points_coordinates_it = interp_points_coord.begin(nb_interpolation_points_per_elem, nb_quad_per_element); Array::matrix_iterator result_begin = result_vec.begin_reinterpret(field_vec.getNbComponent(), nb_interpolation_points_per_elem, result_vec.getSize() / nb_interpolation_points_per_elem); Array::const_matrix_iterator inv_quad_coord_it = interp_inv_coord.begin(nb_quad_per_element, nb_quad_per_element); /// loop over the elements of the current material and element type for (UInt el = 0; el < nb_element; ++el, ++field_it, ++inv_quad_coord_it, ++interpolation_points_coordinates_it) { /** * matrix containing the inversion of the quadrature points' * coordinates */ const Matrix & inv_quad_coord_matrix = *inv_quad_coord_it; /** * multiply it by the field values over quadrature points to get * the interpolation coefficients */ coefficients.mul(inv_quad_coord_matrix, *field_it); /// matrix containing the points' coordinates const Matrix & coord = *interpolation_points_coordinates_it; /// multiply the coordinates matrix by the coefficients matrix and store the result Matrix res(result_begin[elem_fil(el)]); res.mul(coefficients, coord); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::interpolateStressOnFacets(ElementTypeMapArray & result, const GhostType ghost_type) { interpolateElementalFieldOnFacets(stress, result, ghost_type); } /* -------------------------------------------------------------------------- */ void Material::interpolateElementalFieldOnFacets(const ElementTypeMapArray & field, ElementTypeMapArray & result, const GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt sp2 = spatial_dimension * spatial_dimension; const Mesh & mesh = this->model->getMesh(); const Mesh & mesh_facets = mesh.getMeshFacets(); Mesh::type_iterator it = element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last = element_filter.lastType(spatial_dimension, ghost_type); for (; it != last; ++it) { ElementType type = *it; Array & elem_fil = element_filter(type, ghost_type); UInt nb_element = elem_fil.getSize(); UInt nb_quad_per_element = fem->getNbQuadraturePoints(type, ghost_type); const Array & field_vec = field(type, ghost_type); Matrix coefficients(nb_quad_per_element, field_vec.getNbComponent()); const Array & interp_inv_coord = interpolation_inverse_coordinates(type, ghost_type); const Array & interp_points_coord = interpolation_points_matrices(type, ghost_type); UInt nb_interpolation_points_per_elem = interp_points_coord.getNbComponent() / nb_quad_per_element; Array::const_matrix_iterator field_it = field_vec.begin_reinterpret(field_vec.getNbComponent(), nb_quad_per_element, nb_element); Array::const_matrix_iterator interpolation_points_coordinates_it = interp_points_coord.begin(nb_interpolation_points_per_elem, nb_quad_per_element); Array::const_matrix_iterator inv_quad_coord_it = interp_inv_coord.begin(nb_quad_per_element, nb_quad_per_element); Matrix result_tmp(sp2, nb_interpolation_points_per_elem); const Array & facet_to_element = mesh_facets.getSubelementToElement(type, ghost_type); ElementType type_facet = Mesh::getFacetType(type); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); UInt nb_quad_per_facet = nb_interpolation_points_per_elem / nb_facet_per_elem; Element element_for_comparison(type, 0, ghost_type); const Array< std::vector > * element_to_facet = NULL; GhostType current_ghost_type = _casper; Array * result_vec = NULL; /// loop over the elements of the current material and element type for (UInt el = 0; el < nb_element; ++el, ++field_it, ++inv_quad_coord_it, ++interpolation_points_coordinates_it) { /** * matrix containing the inversion of the quadrature points' * coordinates */ const Matrix & inv_quad_coord_matrix = *inv_quad_coord_it; /** * multiply it by the field values over quadrature points to get * the interpolation coefficients */ coefficients.mul(inv_quad_coord_matrix, *field_it); /// matrix containing the points' coordinates const Matrix & coord = *interpolation_points_coordinates_it; /// multiply the coordinates matrix by the coefficients matrix and store the result result_tmp.mul(coefficients, coord); UInt global_el = elem_fil(el); element_for_comparison.element = global_el; for (UInt f = 0; f < nb_facet_per_elem; ++f) { Element facet_elem = facet_to_element(global_el, f); UInt global_facet = facet_elem.element; if (facet_elem.ghost_type != current_ghost_type) { current_ghost_type = facet_elem.ghost_type; element_to_facet = &mesh_facets.getElementToSubelement(type_facet, current_ghost_type); result_vec = &result(type_facet, current_ghost_type); } bool is_second_element = (*element_to_facet)(global_facet)[0] != element_for_comparison; for (UInt q = 0; q < nb_quad_per_facet; ++q) { Vector result_local(result_vec->storage() + (global_facet * nb_quad_per_facet + q) * result_vec->getNbComponent() + is_second_element * sp2, sp2); result_local = result_tmp(f * nb_quad_per_facet + q); } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ +template +const Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) const { + AKANTU_DEBUG_TO_IMPLEMENT(); + return NULL; +} +/* -------------------------------------------------------------------------- */ +template +Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) { + AKANTU_DEBUG_TO_IMPLEMENT(); + return NULL; +} + +/* -------------------------------------------------------------------------- */ +template<> const Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) const { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch(debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" < Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch(debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ +template +const InternalField & Material::getInternal(const ID & int_id) const { + AKANTU_DEBUG_TO_IMPLEMENT(); + return NULL; +} + +/* -------------------------------------------------------------------------- */ +template +InternalField & Material::getInternal(const ID & int_id) { + AKANTU_DEBUG_TO_IMPLEMENT(); + return NULL; +} + +/* -------------------------------------------------------------------------- */ +template<> const InternalField & Material::getInternal(const ID & int_id) const { std::map *>::const_iterator it = internal_vectors_real.find(getID() + ":" + int_id); if(it == internal_vectors_real.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ +template<> InternalField & Material::getInternal(const ID & int_id) { std::map *>::iterator it = internal_vectors_real.find(getID() + ":" + int_id); if(it == internal_vectors_real.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ -const InternalField & Material::getInternalUInt(const ID & int_id) const { +template<> +const InternalField & Material::getInternal(const ID & int_id) const { std::map *>::const_iterator it = internal_vectors_uint.find(getID() + ":" + int_id); if(it == internal_vectors_uint.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ -InternalField & Material::getInternalUInt(const ID & int_id) { +template<> +InternalField & Material::getInternal(const ID & int_id) { std::map *>::iterator it = internal_vectors_uint.find(getID() + ":" + int_id); if(it == internal_vectors_uint.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ void Material::addElements(const Array & elements_to_add) { AKANTU_DEBUG_IN(); UInt mat_id = model->getInternalIndexFromID(getID()); Array::const_iterator el_begin = elements_to_add.begin(); Array::const_iterator el_end = elements_to_add.end(); for(;el_begin != el_end; ++el_begin) { const Element & element = *el_begin; Array & mat_indexes = model->getMaterialByElement (element.type, element.ghost_type); Array & mat_loc_num = model->getMaterialLocalNumbering(element.type, element.ghost_type); UInt index = this->addElement(element.type, element.element, element.ghost_type); mat_indexes(element.element) = mat_id; mat_loc_num(element.element) = index; } this->resizeInternals(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::removeElements(const Array & elements_to_remove) { AKANTU_DEBUG_IN(); Array::const_iterator el_begin = elements_to_remove.begin(); Array::const_iterator el_end = elements_to_remove.end(); if(el_begin==el_end) return; ElementTypeMapArray material_local_new_numbering("remove mat filter elem", getID()); Element element; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; element.ghost_type = ghost_type; ElementTypeMapArray::type_iterator it = element_filter.firstType(_all_dimensions, ghost_type, _ek_not_defined); ElementTypeMapArray::type_iterator end = element_filter.lastType(_all_dimensions, ghost_type, _ek_not_defined); for(; it != end; ++it) { ElementType type = *it; element.type = type; Array & elem_filter = this->element_filter(type, ghost_type); Array & mat_loc_num = this->model->getMaterialLocalNumbering(type, ghost_type); if(!material_local_new_numbering.exists(type, ghost_type)) material_local_new_numbering.alloc(elem_filter.getSize(), 1, type, ghost_type); Array & mat_renumbering = material_local_new_numbering(type, ghost_type); UInt nb_element = elem_filter.getSize(); element.kind=(*el_begin).kind; Array elem_filter_tmp; UInt new_id = 0; for (UInt el = 0; el < nb_element; ++el) { element.element = elem_filter(el); if(std::find(el_begin, el_end, element) == el_end) { elem_filter_tmp.push_back(element.element); mat_renumbering(el) = new_id; mat_loc_num(element.element) = new_id; ++new_id; } else { mat_renumbering(el) = UInt(-1); } } elem_filter.resize(elem_filter_tmp.getSize()); elem_filter.copy(elem_filter_tmp); } } for (std::map *>::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->removeQuadraturePoints(material_local_new_numbering); for (std::map *>::iterator it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->removeQuadraturePoints(material_local_new_numbering); for (std::map *>::iterator it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->removeQuadraturePoints(material_local_new_numbering); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::resizeInternals() { AKANTU_DEBUG_IN(); for (std::map *>::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->resize(); for (std::map *>::iterator it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->resize(); for (std::map *>::iterator it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->resize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::onElementsAdded(__attribute__((unused)) const Array & element_list, __attribute__((unused)) const NewElementsEvent & event) { this->resizeInternals(); } /* -------------------------------------------------------------------------- */ void Material::onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, __attribute__((unused)) const RemovedElementsEvent & event) { UInt my_num = model->getInternalIndexFromID(getID()); ElementTypeMapArray material_local_new_numbering("remove mat filter elem", getID()); Array::const_iterator el_begin = element_list.begin(); Array::const_iterator el_end = element_list.end(); for (ghost_type_t::iterator g = ghost_type_t::begin(); g != ghost_type_t::end(); ++g) { GhostType gt = *g; ElementTypeMapArray::type_iterator it = new_numbering.firstType(_all_dimensions, gt, _ek_not_defined); ElementTypeMapArray::type_iterator end = new_numbering.lastType (_all_dimensions, gt, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; if(element_filter.exists(type, gt) && element_filter(type, gt).getSize()){ Array & elem_filter = element_filter(type, gt); Array & mat_indexes = this->model->getMaterialByElement (*it, gt); Array & mat_loc_num = this->model->getMaterialLocalNumbering(*it, gt); UInt nb_element = this->model->getMesh().getNbElement(type, gt); // all materials will resize of the same size... mat_indexes.resize(nb_element); mat_loc_num.resize(nb_element); if(!material_local_new_numbering.exists(type, gt)) material_local_new_numbering.alloc(elem_filter.getSize(), 1, type, gt); Array & mat_renumbering = material_local_new_numbering(type, gt); const Array & renumbering = new_numbering(type, gt); Array elem_filter_tmp; UInt ni = 0; Element el; el.type = type; el.ghost_type = gt; + el.kind = Mesh::getKind(type); for (UInt i = 0; i < elem_filter.getSize(); ++i) { el.element = elem_filter(i); if(std::find(el_begin, el_end, el) == el_end) { UInt new_el = renumbering(el.element); AKANTU_DEBUG_ASSERT(new_el != UInt(-1), "A not removed element as been badly renumbered"); elem_filter_tmp.push_back(new_el); mat_renumbering(i) = ni; mat_indexes(new_el) = my_num; mat_loc_num(new_el) = ni; ++ni; } else { mat_renumbering(i) = UInt(-1); } } elem_filter.resize(elem_filter_tmp.getSize()); elem_filter.copy(elem_filter_tmp); } } } for (std::map *>::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->removeQuadraturePoints(material_local_new_numbering); for (std::map *>::iterator it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->removeQuadraturePoints(material_local_new_numbering); for (std::map *>::iterator it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->removeQuadraturePoints(material_local_new_numbering); } /* -------------------------------------------------------------------------- */ void Material::onBeginningSolveStep(const AnalysisMethod & method) { this->savePreviousState(); } /* -------------------------------------------------------------------------- */ void Material::onEndSolveStep(const AnalysisMethod & method) { ElementTypeMapArray::type_iterator it = this->element_filter.firstType(_all_dimensions, _not_ghost, _ek_not_defined); ElementTypeMapArray::type_iterator end = element_filter.lastType(_all_dimensions, _not_ghost, _ek_not_defined); for(; it != end; ++it) { this->updateEnergies(*it, _not_ghost); } } /* -------------------------------------------------------------------------- */ void Material::onDamageIteration() { this->savePreviousState(); } /* -------------------------------------------------------------------------- */ void Material::onDamageUpdate() { ElementTypeMapArray::type_iterator it = this->element_filter.firstType(_all_dimensions, _not_ghost, _ek_not_defined); ElementTypeMapArray::type_iterator end = element_filter.lastType(_all_dimensions, _not_ghost, _ek_not_defined); for(; it != end; ++it) { if(!this->potential_energy.exists(*it, _not_ghost)) { UInt nb_element = this->element_filter(*it, _not_ghost).getSize(); UInt nb_quadrature_points = this->fem->getNbQuadraturePoints(*it, _not_ghost); this->potential_energy.alloc(nb_element * nb_quadrature_points, 1, *it, _not_ghost); } this->updateEnergiesAfterDamage(*it, _not_ghost); } } /* -------------------------------------------------------------------------- */ void Material::onDump(){ if(this->isFiniteDeformation()) this->computeAllCauchyStresses(_not_ghost); } /* -------------------------------------------------------------------------- */ void Material::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); std::string type = getID().substr(getID().find_last_of(":") + 1); stream << space << "Material " << type << " [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ inline ElementTypeMap Material::getInternalDataPerElem(const ID & id, const ElementKind & element_kind, const ID & fe_engine_id) const { std::map *>::const_iterator internal_array = internal_vectors_real.find(this->getID()+":"+id); if (internal_array == internal_vectors_real.end() || internal_array->second->getElementKind() != element_kind) AKANTU_EXCEPTION("Cannot find internal field " << id << " in material " << this->name); InternalField & internal = *internal_array->second; InternalField::type_iterator it = internal.firstType(internal.getSpatialDimension(), _not_ghost, element_kind); InternalField::type_iterator last_type = internal.lastType(internal.getSpatialDimension(), _not_ghost, element_kind); ElementTypeMap res; for(; it != last_type; ++it) { UInt nb_quadrature_points = 0; nb_quadrature_points = model->getFEEngine(fe_engine_id).getNbQuadraturePoints(*it); res(*it) = internal.getNbComponent() * nb_quadrature_points; } return res; } /* -------------------------------------------------------------------------- */ void Material::flattenInternal(const std::string & field_id, ElementTypeMapArray & internal_flat, const GhostType ghost_type, ElementKind element_kind) const { this->flattenInternalIntern(field_id, internal_flat, this->spatial_dimension, ghost_type, element_kind); } /* -------------------------------------------------------------------------- */ void Material::flattenInternalIntern(const std::string & field_id, ElementTypeMapArray & internal_flat, UInt spatial_dimension, const GhostType ghost_type, ElementKind element_kind, const ElementTypeMapArray * element_filter, const Mesh * mesh) const { typedef ElementTypeMapArray::type_iterator iterator; if(element_filter == NULL) element_filter = &(this->element_filter); if(mesh == NULL) mesh = &(this->model->mesh); iterator tit = element_filter->firstType(spatial_dimension, ghost_type, element_kind); iterator end = element_filter->lastType(spatial_dimension, ghost_type, element_kind); for (; tit != end; ++tit) { ElementType type = *tit; try { __attribute__((unused)) const Array & src_vect - = this->getArray(field_id, type, ghost_type); + = this->getArray(field_id, type, ghost_type); } catch(debug::Exception & e) { continue; } - const Array & src_vect = this->getArray(field_id, type, ghost_type); + const Array & src_vect = this->getArray(field_id, type, ghost_type); const Array & filter = (*element_filter)(type, ghost_type); // total number of elements for a given type UInt nb_element = mesh->getNbElement(type,ghost_type); // number of filtered elements UInt nb_element_src = filter.getSize(); // number of quadrature points per elem UInt nb_quad_per_elem = 0; // number of data per quadrature point UInt nb_data_per_quad = src_vect.getNbComponent(); if (!internal_flat.exists(type,ghost_type)) { internal_flat.alloc(nb_element * nb_quad_per_elem, nb_data_per_quad, type, ghost_type); } if (nb_element_src == 0) continue; nb_quad_per_elem = (src_vect.getSize() / nb_element_src); // number of data per element UInt nb_data = nb_quad_per_elem * src_vect.getNbComponent(); Array & dst_vect = internal_flat(type,ghost_type); dst_vect.resize(nb_element*nb_quad_per_elem); Array::const_scalar_iterator it = filter.begin(); Array::const_scalar_iterator end = filter.end(); Array::const_vector_iterator it_src = src_vect.begin_reinterpret(nb_data, nb_element_src); Array::vector_iterator it_dst = dst_vect.begin_reinterpret(nb_data, nb_element); for (; it != end ; ++it,++it_src) { it_dst[*it] = *it_src; } } }; /* -------------------------------------------------------------------------- */ - +/// extrapolate internal values +void Material::extrapolateInternal(const ID & id, const Element & element, const Matrix & point, Matrix & extrapolated) { + if (this->isInternal(id, element.kind)) { + UInt nb_element = this->element_filter(element.type, element.ghost_type).getSize(); + const ID name = this->getID() + ":" + id; + UInt nb_quads = this->internal_vectors_real[name]->getFEEngine().getNbQuadraturePoints(element.type, element.ghost_type); + const Array & internal = this->getArray(id, element.type, element.ghost_type); + UInt nb_component = internal.getNbComponent(); + Array::const_matrix_iterator internal_it = internal.begin_reinterpret(nb_component, nb_quads, nb_element); + Element local_element = this->convertToLocalElement(element); + + /// instead of really extrapolating, here the value of the first GP + /// is copied into the result vector. This works only for linear + /// elements + /// @todo extrapolate!!!! + + const Matrix & values = internal_it[local_element.element]; + UInt index = 0; + Vector tmp(nb_component); + for (UInt j = 0; j < values.cols(); ++j) { + tmp = values(j); + if (tmp.norm() > 0) { + index = j; + continue; + } + } + + for (UInt i = 0; i < extrapolated.size(); ++i) { + extrapolated(i) = values(index); + } + } + else { + Matrix default_values(extrapolated.rows(), extrapolated.cols(), 0.); + extrapolated = default_values; + } +} __END_AKANTU__ diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh index 28acc2817..a5634a5f8 100644 --- a/src/model/solid_mechanics/material.hh +++ b/src/model/solid_mechanics/material.hh @@ -1,638 +1,645 @@ /** * @file material.hh * * @author Marco Vocialta * @author Nicolas Richart * @author Daniel Pino Muñoz * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Sep 16 2014 * * @brief Mother class for all materials * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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. * * Akantu 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "aka_voigthelper.hh" #include "parser.hh" #include "parsable.hh" #include "data_accessor.hh" #include "internal_field.hh" #include "random_internal_field.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_HH__ #define __AKANTU_MATERIAL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class SolidMechanicsModel; } __BEGIN_AKANTU__ /** * Interface of all materials * Prerequisites for a new material * - inherit from this class * - implement the following methods: * \code * virtual Real getStableTimeStep(Real h, const Element & element = ElementNull); * * virtual void computeStress(ElementType el_type, * GhostType ghost_type = _not_ghost); * * virtual void computeTangentStiffness(const ElementType & el_type, * Array & tangent_matrix, * GhostType ghost_type = _not_ghost); * \endcode * */ class Material : public Memory, public DataAccessor, public Parsable, public MeshEventHandler, protected SolidMechanicsModelEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /// Initialize material with defaults Material(SolidMechanicsModel & model, const ID & id = ""); /// Initialize material with custom mesh & fe_engine Material(SolidMechanicsModel & model, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id = ""); /// Destructor virtual ~Material(); protected: void initialize(); /* ------------------------------------------------------------------------ */ /* Function that materials can/should reimplement */ /* ------------------------------------------------------------------------ */ protected: /// constitutive law virtual void computeStress(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the tangent stiffness matrix virtual void computeTangentModuli(__attribute__((unused)) const ElementType & el_type, __attribute__((unused)) Array & tangent_matrix, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the potential energy virtual void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the potential energy for an element virtual void computePotentialEnergyByElement(__attribute__((unused)) ElementType type, __attribute__((unused)) UInt index, __attribute__((unused)) Vector & epot_on_quad_points) { AKANTU_DEBUG_TO_IMPLEMENT(); } virtual void updateEnergies(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) { } virtual void updateEnergiesAfterDamage(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) {} /// set the material to steady state (to be implemented for materials that need it) virtual void setToSteadyState(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) { } /// function called to update the internal parameters when the modifiable /// parameters are modified virtual void updateInternalParameters() {} public: + /// extrapolate internal values + virtual void extrapolateInternal(const ID & id, const Element & element, const Matrix & points, Matrix & extrapolated); /// compute the p-wave speed in the material virtual Real getPushWaveSpeed(const Element & element) const { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the s-wave speed in the material virtual Real getShearWaveSpeed(const Element & element) const { AKANTU_DEBUG_TO_IMPLEMENT(); } /// get a material celerity to compute the stable time step (default: is the push wave speed) virtual Real getCelerity(const Element & element) const { return getPushWaveSpeed(element); } /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: template void registerInternal(__attribute__((unused)) InternalField & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); } template void unregisterInternal(__attribute__((unused)) InternalField & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// initialize the material computed parameter virtual void initMaterial(); /// compute the residual for this material virtual void updateResidual(GhostType ghost_type = _not_ghost); /// assemble the residual for this material virtual void assembleResidual(GhostType ghost_type); /// Operations before and after solveStep in implicit virtual void beforeSolveStep() {} virtual void afterSolveStep() {} /// save the stress in the previous_stress if needed virtual void savePreviousState(); /// compute the stresses for this material virtual void computeAllStresses(GhostType ghost_type = _not_ghost); virtual void computeAllNonLocalStresses(__attribute__((unused)) GhostType ghost_type = _not_ghost) {}; virtual void computeAllStressesFromTangentModuli(GhostType ghost_type = _not_ghost); virtual void computeAllCauchyStresses(GhostType ghost_type = _not_ghost); /// set material to steady state void setToSteadyState(GhostType ghost_type = _not_ghost); /// compute the stiffness matrix virtual void assembleStiffnessMatrix(GhostType ghost_type); /// add an element to the local mesh filter inline UInt addElement(const ElementType & type, UInt element, const GhostType & ghost_type); /// add many elements at once void addElements(const Array & elements_to_add); /// remove many element at once void removeElements(const Array & elements_to_remove); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /** * interpolate stress on given positions for each element by means * of a geometrical interpolation on quadrature points */ void interpolateStress(ElementTypeMapArray & result, const GhostType ghost_type = _not_ghost); /** * interpolate stress on given positions for each element by means * of a geometrical interpolation on quadrature points and store the * results per facet */ void interpolateStressOnFacets(ElementTypeMapArray & result, const GhostType ghost_type = _not_ghost); /** * function to initialize the elemental field interpolation * function by inverting the quadrature points' coordinates */ void initElementalFieldInterpolation(const ElementTypeMapArray & interpolation_points_coordinates); /* ------------------------------------------------------------------------ */ /* Common part */ /* ------------------------------------------------------------------------ */ protected: /// assemble the residual template void assembleResidual(GhostType ghost_type); /// Computation of Cauchy stress tensor in the case of finite deformation template void computeCauchyStress(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost); template inline void computeCauchyStressOnQuad(const Matrix & F, const Matrix & S, Matrix & cauchy, const Real & C33 = 1.0 ) const; template void computeAllStressesFromTangentModuli(const ElementType & type, GhostType ghost_type); template void assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type); /// assembling in finite deformation template void assembleStiffnessMatrixNL(const ElementType & type, GhostType ghost_type); template void assembleStiffnessMatrixL2(const ElementType & type, GhostType ghost_type); /// write the stress tensor in the Voigt notation. template inline void SetCauchyStressArray(const Matrix & S_t, Matrix & Stress_vect); inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const; /// Size of the Stress matrix for the case of finite deformation see: Bathe et al, IJNME, Vol 9, 353-386, 1975 inline UInt getCauchyStressMatrixSize(UInt spatial_dimension) const; /// Sets the stress matrix according to Bathe et al, IJNME, Vol 9, 353-386, 1975 template inline void setCauchyStressMatrix(const Matrix & S_t, Matrix & Stress_matrix); /// compute the potential energy by element void computePotentialEnergyByElements(); /// resize the intenals arrays virtual void resizeInternals(); public: /// compute the coordinates of the quadrature points virtual void computeQuadraturePointsCoordinates(ElementTypeMapArray & quadrature_points_coordinates, const GhostType & ghost_type) const; protected: /// interpolate an elemental field on given points for each element void interpolateElementalField(const ElementTypeMapArray & field, ElementTypeMapArray & result, const GhostType ghost_type); /// interpolate an elemental field and store the results per facet void interpolateElementalFieldOnFacets(const ElementTypeMapArray & field, ElementTypeMapArray & result, const GhostType ghost_type); /// template function to initialize the elemental field interpolation template void initElementalFieldInterpolation(const Array & quad_coordinates, const Array & interpolation_points_coordinates, const UInt nb_interpolation_points_per_elem, const GhostType ghost_type); /// build the coordinate matrix for the interpolation on elemental field template inline void buildElementalFieldInterpolationCoodinates(const Matrix & coordinates, Matrix & coordMatrix); /// build interpolation coordinates for basic linear elements inline void buildElementalFieldInterpolationCoodinatesLinear(const Matrix & coordinates, Matrix & coordMatrix); /// build interpolation coordinates for basic quadratic elements inline void buildElementalFieldInterpolationCoodinatesQuadratic(const Matrix & coordinates, Matrix & coordMatrix); public: /* ------------------------------------------------------------------------ */ /* Conversion functions */ /* ------------------------------------------------------------------------ */ template static inline void gradUToF (const Matrix & grad_u, Matrix & F); static inline void rightCauchy(const Matrix & F, Matrix & C); static inline void leftCauchy (const Matrix & F, Matrix & B); template static inline void gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon); template static inline void gradUToGreenStrain(const Matrix & grad_u, Matrix & epsilon); static inline Real stressToVonMises(const Matrix & stress); protected: /// converts global element to local element inline Element convertToLocalElement(const Element & global_element) const; /// converts local element to global element inline Element convertToGlobalElement(const Element & local_element) const; /// converts global quadrature point to local quadrature point inline QuadraturePoint convertToLocalPoint(const QuadraturePoint & global_point) const; /// converts local quadrature point to global quadrature point inline QuadraturePoint convertToGlobalPoint(const QuadraturePoint & local_point) const; /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: virtual inline UInt getNbDataForElements(const Array & elements, SynchronizationTag tag) const; virtual inline void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const; virtual inline void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag); template inline void packElementDataHelper(const ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()) const; template inline void unpackElementDataHelper(ElementTypeMapArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()); /* ------------------------------------------------------------------------ */ /* MeshEventHandler inherited members */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ virtual void onElementsAdded(const Array & element_list, const NewElementsEvent & event); virtual void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event); /* ------------------------------------------------------------------------ */ /* SolidMechanicsModelEventHandler inherited members */ /* ------------------------------------------------------------------------ */ public: virtual void onBeginningSolveStep(const AnalysisMethod & method); virtual void onEndSolveStep(const AnalysisMethod & method); virtual void onDamageIteration(); virtual void onDamageUpdate(); virtual void onDump(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Name, name, const std::string &); AKANTU_GET_MACRO(Model, *model, const SolidMechanicsModel &) AKANTU_GET_MACRO(ID, Memory::getID(), const ID &); AKANTU_GET_MACRO(Rho, rho, Real); AKANTU_SET_MACRO(Rho, rho, Real); AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// return the potential energy for the subset of elements contained by the material Real getPotentialEnergy(); /// return the potential energy for the provided element Real getPotentialEnergy(ElementType & type, UInt index); /// return the energy (identified by id) for the subset of elements contained by the material virtual Real getEnergy(std::string energy_id); /// return the energy (identified by id) for the provided element virtual Real getEnergy(std::string energy_id, ElementType type, UInt index); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(GradU, gradu, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy, Real); AKANTU_GET_MACRO(GradU, gradu, const ElementTypeMapArray &); AKANTU_GET_MACRO(Stress, stress, const ElementTypeMapArray &); AKANTU_GET_MACRO(ElementFilter, element_filter, const ElementTypeMapArray &); AKANTU_GET_MACRO(FEEngine, *fem, FEEngine &); bool isNonLocal() const { return is_non_local; } - const Array & getArray(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost) const; - Array & getArray(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost); - - const InternalField & getInternal(const ID & id) const; - InternalField & getInternal(const ID & id); + template + const Array & getArray(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost) const; + template + Array & getArray(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost); - const InternalField & getInternalUInt(const ID & id) const; - InternalField & getInternalUInt(const ID & id); + template + const InternalField & getInternal(const ID & id) const; + template + InternalField & getInternal(const ID & id); inline bool isInternal(const ID & id, const ElementKind & element_kind) const; virtual ElementTypeMap getInternalDataPerElem(const ID & id, const ElementKind & element_kind, const ID & fe_engine_id = "") const; bool isFiniteDeformation() const { return finite_deformation; } bool isInelasticDeformation() const { return inelastic_deformation; } template inline void setParam(const ID & param, T value); template inline const T & getParam(const ID & param) const; virtual void flattenInternal(const std::string & field_id, ElementTypeMapArray & internal_flat, const GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined) const; protected: /// internal variation of the flatten function that is more flexible and can /// be used by inherited materials to change some behavior virtual void flattenInternalIntern(const std::string & field_id, ElementTypeMapArray & internal_flat, UInt spatial_dimension, const GhostType ghost_type, ElementKind element_kind, const ElementTypeMapArray * element_filter = NULL, const Mesh * mesh = NULL) const; protected: bool isInit() const { return is_init; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ -private: +protected: /// boolean to know if the material has been initialized bool is_init; std::map *> internal_vectors_real; std::map *> internal_vectors_uint; std::map *> internal_vectors_bool; protected: /// Link to the fem object in the model FEEngine * fem; /// Finite deformation bool finite_deformation; /// Finite deformation bool inelastic_deformation; /// material name std::string name; /// The model to witch the material belong SolidMechanicsModel * model; /// density : rho Real rho; /// spatial dimension UInt spatial_dimension; /// list of element handled by the material ElementTypeMapArray element_filter; /// stresses arrays ordered by element types InternalField stress; /// eigenstrain arrays ordered by element types InternalField eigenstrain; /// grad_u arrays ordered by element types InternalField gradu; /// Green Lagrange strain (Finite deformation) InternalField green_strain; /// Second Piola-Kirchhoff stress tensor arrays ordered by element types (Finite deformation) InternalField piola_kirchhoff_2; /// potential energy by element InternalField potential_energy; /// tell if using in non local mode or not bool is_non_local; /// tell if the material need the previous stress state bool use_previous_stress; /// tell if the material need the previous strain state bool use_previous_gradu; /// elemental field interpolation coordinates InternalField interpolation_inverse_coordinates; /// elemental field interpolation points InternalField interpolation_points_matrices; + + /// vector that contains the names of all the internals that need to + /// be transferred when material interfaces move + std::vector internals_to_transfer; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Material & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #include "internal_field_tmpl.hh" #include "random_internal_field_tmpl.hh" /* -------------------------------------------------------------------------- */ /* Auto loop */ /* -------------------------------------------------------------------------- */ #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type) \ Array::matrix_iterator gradu_it = \ this->gradu(el_type, ghost_type).begin(this->spatial_dimension, \ this->spatial_dimension); \ Array::matrix_iterator gradu_end = \ this->gradu(el_type, ghost_type).end(this->spatial_dimension, \ this->spatial_dimension); \ \ this->stress(el_type, \ ghost_type).resize(this->gradu(el_type, \ ghost_type).getSize()); \ \ Array::iterator< Matrix > stress_it = \ this->stress(el_type, ghost_type).begin(this->spatial_dimension, \ this->spatial_dimension); \ \ if(this->isFiniteDeformation()){ \ this->piola_kirchhoff_2(el_type, \ ghost_type).resize(this->gradu(el_type, \ ghost_type).getSize()); \ stress_it = \ this->piola_kirchhoff_2(el_type, \ ghost_type).begin(this->spatial_dimension, \ this->spatial_dimension); \ } \ \ for(;gradu_it != gradu_end; ++gradu_it, ++stress_it) { \ Matrix & __attribute__((unused)) grad_u = *gradu_it; \ Matrix & __attribute__((unused)) sigma = *stress_it #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END \ } \ #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \ Array::matrix_iterator gradu_it = \ this->gradu(el_type, ghost_type).begin(this->spatial_dimension, \ this->spatial_dimension); \ Array::matrix_iterator gradu_end = \ this->gradu(el_type, ghost_type).end(this->spatial_dimension, \ this->spatial_dimension); \ Array::matrix_iterator sigma_it = \ this->stress(el_type, ghost_type).begin(this->spatial_dimension, \ this->spatial_dimension); \ \ tangent_mat.resize(this->gradu(el_type, ghost_type).getSize()); \ \ UInt tangent_size = \ this->getTangentStiffnessVoigtSize(this->spatial_dimension); \ Array::matrix_iterator tangent_it = \ tangent_mat.begin(tangent_size, \ tangent_size); \ \ for(;gradu_it != gradu_end; ++gradu_it, ++sigma_it, ++tangent_it) { \ Matrix & __attribute__((unused)) grad_u = *gradu_it; \ Matrix & __attribute__((unused)) sigma_tensor = *sigma_it; \ Matrix & tangent = *tangent_it #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END \ } \ /* -------------------------------------------------------------------------- */ #define INSTANTIATE_MATERIAL(mat_name) \ template class mat_name<1>; \ template class mat_name<2>; \ template class mat_name<3> #endif /* __AKANTU_MATERIAL_HH__ */ diff --git a/src/model/solid_mechanics/materials/internal_field.hh b/src/model/solid_mechanics/materials/internal_field.hh index 83a99b092..52de6e0ca 100644 --- a/src/model/solid_mechanics/materials/internal_field.hh +++ b/src/model/solid_mechanics/materials/internal_field.hh @@ -1,215 +1,217 @@ /** * @file internal_field.hh * * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Tue Sep 02 2014 * * @brief Material internal properties * * @section LICENSE * * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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. * * Akantu 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "element_type_map.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_INTERNAL_FIELD_HH__ #define __AKANTU_INTERNAL_FIELD_HH__ __BEGIN_AKANTU__ class Material; class FEEngine; template class InternalField : public ElementTypeMapArray { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: InternalField(const ID & id, Material & material); virtual ~InternalField(); /// This constructor is only here to let cohesive elements compile InternalField(const ID & id, Material & material, FEEngine & fem, const ElementTypeMapArray & element_filter); /// More general constructor InternalField(const ID & id, Material & material, UInt dim, FEEngine & fem, const ElementTypeMapArray & element_filter); InternalField(const ID & id, const InternalField & other); private: InternalField operator=(__attribute__((unused)) const InternalField & other) {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// function to reset the FEEngine for the internal field virtual void setFEEngine(FEEngine & fe_engine); /// function to reset the element kind for the internal virtual void setElementKind(ElementKind element_kind); /// initialize the field to a given number of component virtual void initialize(UInt nb_component); /// activate the history of this field virtual void initializeHistory(); /// resize the arrays and set the new element to 0 virtual void resize(); /// set the field to a given value v virtual void setDefaultValue(const T & v); /// reset all the fields to the default value virtual void reset(); /// save the current values in the history virtual void saveCurrentValues(); /// remove the quadrature points corresponding to suppressed elements virtual void removeQuadraturePoints(const ElementTypeMapArray & new_numbering); /// print the content virtual void printself(std::ostream & stream, UInt indent = 0) const; /// get the default value inline operator T() const; - AKANTU_GET_MACRO(FEEngine, *fem, FEEngine &); + virtual FEEngine & getFEEngine() { + return *fem;} + ///AKANTU_GET_MACRO(FEEngine, *fem, FEEngine &); protected: /// initialize the arrays in the ElementTypeMapArray void internalInitialize(UInt nb_component); /// set the values for new internals virtual void setArrayValues(T * begin, T * end); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the Array corresponding to the type en ghost_type specified virtual Array & operator()(const ElementType & type, const GhostType & ghost_type = _not_ghost) { return ElementTypeMapArray::operator()(type, ghost_type); } virtual const Array & operator()(const ElementType & type, const GhostType & ghost_type = _not_ghost) const { return ElementTypeMapArray::operator()(type, ghost_type); } virtual Array & previous(const ElementType & type, const GhostType & ghost_type = _not_ghost) { AKANTU_DEBUG_ASSERT(previous_values != NULL, "The history of the internal " << this->getID() << " has not been activated"); return this->previous_values->operator()(type, ghost_type); } virtual const Array & previous(const ElementType & type, const GhostType & ghost_type = _not_ghost) const { AKANTU_DEBUG_ASSERT(previous_values != NULL, "The history of the internal " << this->getID() << " has not been activated"); return this->previous_values->operator()(type, ghost_type); } virtual InternalField & previous() { AKANTU_DEBUG_ASSERT(previous_values != NULL, "The history of the internal " << this->getID() << " has not been activated"); return *(this->previous_values); } virtual const InternalField & previous() const { AKANTU_DEBUG_ASSERT(previous_values != NULL, "The history of the internal " << this->getID() << " has not been activated"); return *(this->previous_values); } /// check if the history is used or not bool hasHistory() const { return (previous_values != NULL); } /// get the kind treated by the internal const ElementKind & getElementKind() const { return element_kind; } /// return the number of components UInt getNbComponent() const { return nb_component; } /// return the spatial dimension corresponding to the internal element type loop filter UInt getSpatialDimension() const { return this->spatial_dimension; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// the material for which this is an internal parameter Material & material; /// the fem containing the mesh and the element informations FEEngine * fem; /// Element filter if needed const ElementTypeMapArray & element_filter; /// default value T default_value; /// spatial dimension of the element to consider UInt spatial_dimension; /// ElementKind of the element to consider ElementKind element_kind; /// Number of component of the internal field UInt nb_component; /// Is the field initialized bool is_init; /// previous values InternalField * previous_values; }; /// standard output stream operator template inline std::ostream & operator <<(std::ostream & stream, const InternalField & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_INTERNAL_FIELD_HH__ */ diff --git a/src/model/solid_mechanics/materials/weight_function.hh b/src/model/solid_mechanics/materials/weight_function.hh index 8a8b5b07c..0400e059d 100644 --- a/src/model/solid_mechanics/materials/weight_function.hh +++ b/src/model/solid_mechanics/materials/weight_function.hh @@ -1,374 +1,374 @@ /** * @file weight_function.hh * * @author Nicolas Richart * @author Cyprien Wolff * * @date creation: Fri Apr 13 2012 * @date last modification: Thu Jun 05 2014 * * @brief Weight functions for non local materials * * @section LICENSE * * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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. * * Akantu 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_types.hh" #include "solid_mechanics_model.hh" #include "parsable.hh" #include #if defined(AKANTU_DEBUG_TOOLS) #include "aka_debug_tools.hh" #include #endif #include "material_list.hh" /* -------------------------------------------------------------------------- */ #include #include "material_damage.hh" #ifndef __AKANTU_WEIGHT_FUNCTION_HH__ #define __AKANTU_WEIGHT_FUNCTION_HH__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Normal weight function */ /* -------------------------------------------------------------------------- */ template class BaseWeightFunction : public Parsable { public: BaseWeightFunction(Material & material, const std::string & type = "base") : Parsable(_st_non_local, "weight_function:" + type), material(material), type(type) { this->registerParam("radius" , R , 100., _pat_parsable | _pat_readable , "Non local radius"); this->registerParam("update_rate" , update_rate, 0U , _pat_parsmod, "Update frequency"); } virtual ~BaseWeightFunction() {} virtual void init() { R2 = R * R; }; virtual void updateInternals(__attribute__((unused)) const ElementTypeMapArray & quadrature_points_coordinates) {}; /* ------------------------------------------------------------------------ */ inline void setRadius(Real radius) { R = radius; R2 = R * R; } /* ------------------------------------------------------------------------ */ inline void selectType(__attribute__((unused)) ElementType type1, __attribute__((unused)) GhostType ghost_type1, __attribute__((unused)) ElementType type2, __attribute__((unused)) GhostType ghost_type2) { } /* ------------------------------------------------------------------------ */ inline Real operator()(Real r, const __attribute__((unused)) QuadraturePoint & q1, const __attribute__((unused)) QuadraturePoint & q2) { Real w = 0; if(r <= R) { Real alpha = (1. - r*r / R2); w = alpha * alpha; // *weight = 1 - sqrt(r / radius); } return w; } void printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "WeightFunction " << type << " [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; } public: Real getRadius() { return R; } UInt getUpdateRate() { return update_rate; } public: virtual UInt getNbDataForElements(__attribute__((unused)) const Array & elements, __attribute__((unused)) SynchronizationTag tag) const { return 0; } virtual inline void packElementData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const Array & elements, __attribute__((unused)) SynchronizationTag tag) const {} virtual inline void unpackElementData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const Array & elements, __attribute__((unused)) SynchronizationTag tag) {} protected: Material & material; Real R; Real R2; UInt update_rate; const std::string type; }; /* -------------------------------------------------------------------------- */ /* Damage weight function */ /* -------------------------------------------------------------------------- */ template class DamagedWeightFunction : public BaseWeightFunction { public: DamagedWeightFunction(Material & material) : BaseWeightFunction(material, "damaged") { //AKANTU_DEBUG_ASSERT(dynamic_cast *>(&material) != NULL, "This weight function works only with damage materials!"); } inline void selectType(__attribute__((unused)) ElementType type1, __attribute__((unused)) GhostType ghost_type1, ElementType type2, GhostType ghost_type2) { - selected_damage = &(this->material.getArray("damage", type2, ghost_type2)); + selected_damage = &(this->material.template getArray("damage", type2, ghost_type2)); } inline Real operator()(Real r, const __attribute__((unused)) QuadraturePoint & q1, const QuadraturePoint & q2) { UInt quad = q2.global_num; Real D = (*selected_damage)(quad); Real Radius_t = 0; Real Radius_init = this->R2; // if(D <= 0.5) // { // Radius_t = 2*D*Radius_init; // } // else // { // Radius_t = 2*Radius_init*(1-D); // } // Radius_t = Radius_init*(1-D); Radius_init *= Radius_init; Radius_t *= Radius_t; if(Radius_t < Math::getTolerance()) { Radius_t = 0.001*Radius_init; } Real expb = (2*std::log(0.51))/(std::log(1.0-0.49*Radius_t/Radius_init)); Int expb_floor=std::floor(expb); Real b = expb_floor + expb_floor%2; Real alpha = std::max(0., 1. - r*r / Radius_init); Real w = std::pow(alpha,b); return w; } private: const Array * selected_damage; }; /* -------------------------------------------------------------------------- */ /* Remove damaged weight function */ /* -------------------------------------------------------------------------- */ template class RemoveDamagedWeightFunction : public BaseWeightFunction { public: RemoveDamagedWeightFunction(Material & material) : BaseWeightFunction(material, "remove_damaged") { this->registerParam("damage_limit", this->damage_limit, 1., _pat_parsable, "Damage Threshold"); } inline void selectType(__attribute__((unused)) ElementType type1, __attribute__((unused)) GhostType ghost_type1, ElementType type2, GhostType ghost_type2) { - selected_damage = &(this->material.getArray("damage", type2, ghost_type2)); + selected_damage = &(this->material.template getArray("damage", type2, ghost_type2)); // selected_damage = &mat.getDamage(type2, ghost_type2); } inline Real operator()(Real r, const __attribute__((unused)) QuadraturePoint & q1, const QuadraturePoint & q2) { UInt quad = q2.global_num; if(q1 == q2) return 1.; Real D = (*selected_damage)(quad); Real w = 0.; if(D < damage_limit) { Real alpha = std::max(0., 1. - r*r / this->R2); w = alpha * alpha; } return w; } virtual UInt getNbDataForElements(const Array & elements, SynchronizationTag tag) const { if(tag == _gst_mnl_weight) return this->material.getModel().getNbQuadraturePoints(elements) * sizeof(Real); return 0; } virtual inline void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const { if(tag == _gst_mnl_weight) { - ElementTypeMapArray & damage = this->material.getInternal("damage"); + ElementTypeMapArray & damage = this->material.template getInternal("damage"); this->material.packElementDataHelper(damage, buffer, elements); } } virtual inline void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) { if(tag == _gst_mnl_weight) { - ElementTypeMapArray & damage = this->material.getInternal("damage"); + ElementTypeMapArray & damage = this->material.template getInternal("damage"); this->material.unpackElementDataHelper(damage, buffer, elements); } } private: /// limit at which a point is considered as complitely broken Real damage_limit; /// internal pointer to the current damage vector const Array * selected_damage; }; /* -------------------------------------------------------------------------- */ /* Remove damaged with damage rate weight function */ /* -------------------------------------------------------------------------- */ template class RemoveDamagedWithDamageRateWeightFunction : public BaseWeightFunction { public: RemoveDamagedWithDamageRateWeightFunction(Material & material) : BaseWeightFunction(material, "remove_damage_with_damage_rate") { this->registerParam("damage_limit", this->damage_limit_with_damage_rate, 1, _pat_parsable, "Damage Threshold"); } inline void selectType(__attribute__((unused)) ElementType type1, __attribute__((unused)) GhostType ghost_type1, ElementType type2, GhostType ghost_type2) { selected_damage_with_damage_rate = &(this->material.getArray("damage",type2, ghost_type2)); selected_damage_rate_with_damage_rate = &(this->material.getArray("damage-rate",type2, ghost_type2)); } inline Real operator()(Real r, const __attribute__((unused)) QuadraturePoint & q1, const QuadraturePoint & q2) { UInt quad = q2.global_num; if(q1.global_num == quad) return 1.; Real D = (*selected_damage_with_damage_rate)(quad); Real w = 0.; Real alphaexp = 1.; Real betaexp = 2.; if(D < damage_limit_with_damage_rate) { Real alpha = std::max(0., 1. - pow((r*r / this->R2),alphaexp)); w = pow(alpha, betaexp); } return w; } private: /// limit at which a point is considered as complitely broken Real damage_limit_with_damage_rate; /// internal pointer to the current damage vector const Array * selected_damage_with_damage_rate; /// internal pointer to the current damage rate vector const Array * selected_damage_rate_with_damage_rate; }; /* -------------------------------------------------------------------------- */ /* Stress Based Weight */ /* -------------------------------------------------------------------------- */ template class StressBasedWeightFunction : public BaseWeightFunction { public: StressBasedWeightFunction(Material & material); void init(); virtual void updateInternals(__attribute__((unused)) const ElementTypeMapArray & quadrature_points_coordinates) { updatePrincipalStress(_not_ghost); updatePrincipalStress(_ghost); }; void updatePrincipalStress(GhostType ghost_type); inline void updateQuadraturePointsCoordinates(ElementTypeMapArray & quadrature_points_coordinates); inline void selectType(ElementType type1, GhostType ghost_type1, ElementType type2, GhostType ghost_type2); inline Real operator()(Real r, const QuadraturePoint & q1, const QuadraturePoint & q2); inline Real computeRhoSquare(Real r, Vector & eigs, Matrix & eigenvects, Vector & x_s); private: Real ft; InternalField stress_diag; Array * selected_stress_diag; InternalField stress_base; Array * selected_stress_base; // InternalField quadrature_points_coordinates; Array * selected_position_1; Array * selected_position_2; InternalField characteristic_size; Array * selected_characteristic_size; }; template inline std::ostream & operator <<(std::ostream & stream, const BaseWeightFunction & _this) { _this.printself(stream); return stream; } #include "weight_function_tmpl.hh" __END_AKANTU__ #endif /* __AKANTU_WEIGHT_FUNCTION_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 6a90777d0..254976275 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1837 +1,1853 @@ /** * @file solid_mechanics_model.cc * * @author Guillaume Anciaux * @author David Simon Kammer * @author Aurelia Isabel Cuba Ramos * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Tue Jul 27 2010 * @date last modification: Fri Sep 19 2014 * * @brief Implementation of the SolidMechanicsModel class * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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. * * Akantu 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_math.hh" #include "aka_common.hh" #include "solid_mechanics_model.hh" #include "group_manager_inline_impl.cc" #include "dumpable_inline_impl.hh" #include "integration_scheme_2nd_order.hh" #include "element_group.hh" #include "static_communicator.hh" #include "dof_synchronizer.hh" #include "element_group.hh" #include #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif #ifdef AKANTU_USE_PETSC #include "solver_petsc.hh" #include "petsc_matrix.hh" #endif #ifdef AKANTU_USE_IOHELPER # include "dumper_field.hh" # include "dumper_paraview.hh" # include "dumper_homogenizing_field.hh" # include "dumper_material_internal_field.hh" # include "dumper_elemental_field.hh" # include "dumper_material_padders.hh" # include "dumper_element_partition.hh" # include "dumper_iohelper.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ const SolidMechanicsModelOptions default_solid_mechanics_model_options(_explicit_lumped_mass, false); /* -------------------------------------------------------------------------- */ /** * A solid mechanics model need a mesh and a dimension to be created. the model * by it self can not do a lot, the good init functions should be called in * order to configure the model depending on what we want to do. * * @param mesh mesh representing the model we want to simulate * @param dim spatial dimension of the problem, if dim = 0 (default value) the * dimension of the problem is assumed to be the on of the mesh * @param id an id to identify the model */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, dim, id, memory_id), BoundaryCondition(), time_step(NAN), f_m2a(1.0), mass_matrix(NULL), velocity_damping_matrix(NULL), stiffness_matrix(NULL), jacobian_matrix(NULL), material_index("material index", id), material_local_numbering("material local numbering", id), material_selector(new DefaultMaterialSelector(material_index)), is_default_material_selector(true), integrator(NULL), increment_flag(false), solver(NULL), synch_parallel(NULL), are_materials_instantiated(false) { AKANTU_DEBUG_IN(); createSynchronizerRegistry(this); registerFEEngineObject("SolidMechanicsFEEngine", mesh, spatial_dimension); this->displacement = NULL; this->mass = NULL; this->velocity = NULL; this->acceleration = NULL; this->force = NULL; this->residual = NULL; this->blocked_dofs = NULL; this->increment = NULL; this->increment_acceleration = NULL; this->dof_synchronizer = NULL; this->previous_displacement = NULL; materials.clear(); mesh.registerEventHandler(*this); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("paraview_all", id, true); this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() { AKANTU_DEBUG_IN(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { delete *mat_it; } materials.clear(); delete integrator; delete solver; delete mass_matrix; delete velocity_damping_matrix; if(stiffness_matrix && stiffness_matrix != jacobian_matrix) delete stiffness_matrix; delete jacobian_matrix; delete synch_parallel; if(is_default_material_selector) { delete material_selector; material_selector = NULL; } AKANTU_DEBUG_OUT(); } void SolidMechanicsModel::setTimeStep(Real time_step) { this->time_step = time_step; #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper().setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ /** * This function groups many of the initialization in on function. For most of * basics case the function should be enough. The functions initialize the * model, the internal vectors, set them to 0, and depending on the parameters * it also initialize the explicit or implicit solver. * * @param material_file the file containing the materials to use * @param method the analysis method wanted. See the akantu::AnalysisMethod for * the different possibilities */ void SolidMechanicsModel::initFull(const ModelOptions & options) { Model::initFull(options); const SolidMechanicsModelOptions & smm_options = dynamic_cast(options); method = smm_options.analysis_method; // initialize the vectors initArrays(); // set the initial condition to 0 force->clear(); velocity->clear(); acceleration->clear(); displacement->clear(); // initialize pcb if(pbc_pair.size()!=0) initPBC(); // initialize the time integration schemes switch(method) { case _explicit_lumped_mass: initExplicit(); break; case _explicit_consistent_mass: initSolver(); initExplicit(); break; case _implicit_dynamic: initImplicit(true); break; case _static: initImplicit(false); break; default: AKANTU_EXCEPTION("analysis method not recognised by SolidMechanicsModel"); break; } // initialize the materials if(this->parser->getLastParsedFile() != "") { instantiateMaterials(); } if(!smm_options.no_init_materials) { initMaterials(); } if(increment_flag) initBC(*this, *displacement, *increment, *force); else initBC(*this, *displacement, *force); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initParallel(MeshPartition * partition, DataAccessor * data_accessor) { AKANTU_DEBUG_IN(); if (data_accessor == NULL) data_accessor = this; synch_parallel = &createParallelSynch(partition,data_accessor); synch_registry->registerSynchronizer(*synch_parallel, _gst_material_id); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_mass); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_stress); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initFEEngineBoundary() { FEEngine & fem_boundary = getFEEngineBoundary(); fem_boundary.initShapeFunctions(_not_ghost); fem_boundary.initShapeFunctions(_ghost); fem_boundary.computeNormalsOnControlPoints(_not_ghost); fem_boundary.computeNormalsOnControlPoints(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initExplicit(AnalysisMethod analysis_method) { AKANTU_DEBUG_IN(); //in case of switch from implicit to explicit if(!this->isExplicit()) method = analysis_method; if (integrator) delete integrator; integrator = new CentralDifference(); UInt nb_nodes = acceleration->getSize(); UInt nb_degree_of_freedom = acceleration->getNbComponent(); std::stringstream sstr; sstr << id << ":increment_acceleration"; increment_acceleration = &(alloc(sstr.str(), nb_nodes, nb_degree_of_freedom, Real())); AKANTU_DEBUG_OUT(); } void SolidMechanicsModel::initArraysPreviousDisplacment() { AKANTU_DEBUG_IN(); SolidMechanicsModel::setIncrementFlagOn(); UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp_t; sstr_disp_t << id << ":previous_displacement"; previous_displacement = &(alloc (sstr_disp_t.str(), nb_nodes, spatial_dimension, 0.)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Allocate all the needed vectors. By default their are not necessarily set to * 0 * */ void SolidMechanicsModel::initArrays() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":displacement"; // std::stringstream sstr_mass; sstr_mass << id << ":mass"; std::stringstream sstr_velo; sstr_velo << id << ":velocity"; std::stringstream sstr_acce; sstr_acce << id << ":acceleration"; std::stringstream sstr_forc; sstr_forc << id << ":force"; std::stringstream sstr_resi; sstr_resi << id << ":residual"; std::stringstream sstr_boun; sstr_boun << id << ":blocked_dofs"; displacement = &(alloc(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); // mass = &(alloc(sstr_mass.str(), nb_nodes, spatial_dimension, 0)); velocity = &(alloc(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); acceleration = &(alloc(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); force = &(alloc(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); residual = &(alloc(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); blocked_dofs = &(alloc(sstr_boun.str(), nb_nodes, spatial_dimension, false)); std::stringstream sstr_curp; sstr_curp << id << ":current_position"; current_position = &(alloc(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE)); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it, gt); material_index.alloc(nb_element, 1, *it, gt); material_local_numbering.alloc(nb_element, 1, *it, gt); } } dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian * */ void SolidMechanicsModel::initModel() { /// \todo add the current position as a parameter to initShapeFunctions for /// large deformation getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initPBC() { Model::initPBC(); registerPBCSynchronizer(); // as long as there are ones on the diagonal of the matrix, we can put boudandary true for slaves std::map::iterator it = pbc_pair.begin(); std::map::iterator end = pbc_pair.end(); UInt dim = mesh.getSpatialDimension(); while(it != end) { for (UInt i=0; iregisterSynchronizer(*synch, _gst_smm_uv); synch_registry->registerSynchronizer(*synch, _gst_smm_mass); synch_registry->registerSynchronizer(*synch, _gst_smm_res); synch_registry->registerSynchronizer(*synch, _gst_for_dump); changeLocalEquationNumberForPBC(pbc_pair, mesh.getSpatialDimension()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); current_position->resize(nb_nodes); Real * current_position_val = current_position->storage(); Real * position_val = mesh.getNodes().storage(); Real * displacement_val = displacement->storage(); /// compute current_position = initial_position + displacement memcpy(current_position_val, position_val, nb_nodes*spatial_dimension*sizeof(Real)); for (UInt n = 0; n < nb_nodes*spatial_dimension; ++n) { *current_position_val++ += *displacement_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initializeUpdateResidualData() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); residual->resize(nb_nodes); /// copy the forces in residual for boundary conditions memcpy(residual->storage(), force->storage(), nb_nodes*spatial_dimension*sizeof(Real)); // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); updateCurrentPosition(); AKANTU_DEBUG_OUT(); } /*----------------------------------------------------------------------------*/ void SolidMechanicsModel::reInitialize() { } /* -------------------------------------------------------------------------- */ /* Explicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * This function compute the second member of the motion equation. That is to * say the sum of forces @f$ r = F_{ext} - F_{int} @f$. @f$ F_{ext} @f$ is given * by the user in the force vector, and @f$ F_{int} @f$ is computed as @f$ * F_{int} = \int_{\Omega} N \sigma d\Omega@f$ * */ void SolidMechanicsModel::updateResidual(bool need_initialize) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); // f = f_ext - f_int // f = f_ext if(need_initialize) initializeUpdateResidualData(); AKANTU_DEBUG_INFO("Compute local stresses"); std::vector::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStresses(_not_ghost); } #ifdef AKANTU_DAMAGE_NON_LOCAL /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ synch_registry->asynchronousSynchronize(_gst_mnl_for_average); AKANTU_DEBUG_INFO("Compute non local stresses for local elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_not_ghost); } AKANTU_DEBUG_INFO("Wait distant non local stresses"); synch_registry->waitEndSynchronize(_gst_mnl_for_average); AKANTU_DEBUG_INFO("Compute non local stresses for ghosts elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_ghost); } #endif /* ------------------------------------------------------------------------ */ /* assembling the forces internal */ // communicate the stress AKANTU_DEBUG_INFO("Send data for residual assembly"); synch_registry->asynchronousSynchronize(_gst_smm_stress); AKANTU_DEBUG_INFO("Assemble residual for local elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.assembleResidual(_not_ghost); } AKANTU_DEBUG_INFO("Wait distant stresses"); // finalize communications synch_registry->waitEndSynchronize(_gst_smm_stress); AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.assembleResidual(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeStresses() { if (isExplicit()) { // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); // compute stresses on all local elements for each materials std::vector::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStresses(_not_ghost); } /* ------------------------------------------------------------------------ */ #ifdef AKANTU_DAMAGE_NON_LOCAL /* Computation of the non local part */ synch_registry->asynchronousSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_not_ghost); } synch_registry->waitEndSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_ghost); } #endif } else { std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStressesFromTangentModuli(_not_ghost); } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateResidualInternal() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Update the residual"); // f = f_ext - f_int - Ma - Cv = r - Ma - Cv; if(method != _static) { // f -= Ma if(mass_matrix) { // if full mass_matrix Array * Ma = new Array(*acceleration, true, "Ma"); *Ma *= *mass_matrix; /// \todo check unit conversion for implicit dynamics // *Ma /= f_m2a *residual -= *Ma; delete Ma; } else if (mass) { // else lumped mass UInt nb_nodes = acceleration->getSize(); UInt nb_degree_of_freedom = acceleration->getNbComponent(); Real * mass_val = mass->storage(); Real * accel_val = acceleration->storage(); Real * res_val = residual->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { if(!(*blocked_dofs_val)) { *res_val -= *accel_val * *mass_val /f_m2a; } blocked_dofs_val++; res_val++; mass_val++; accel_val++; } } else { AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); } // f -= Cv if(velocity_damping_matrix) { Array * Cv = new Array(*velocity); *Cv *= *velocity_damping_matrix; *residual -= *Cv; delete Cv; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateAcceleration() { AKANTU_DEBUG_IN(); updateResidualInternal(); if(method == _explicit_lumped_mass) { /* residual = residual_{n+1} - M * acceleration_n therefore solution = increment acceleration not acceleration */ solveLumped(*increment_acceleration, *mass, *residual, *blocked_dofs, f_m2a); } else if (method == _explicit_consistent_mass) { solve(*increment_acceleration); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveLumped(Array & x, const Array & A, const Array & b, const Array & blocked_dofs, Real alpha) { Real * A_val = A.storage(); Real * b_val = b.storage(); Real * x_val = x.storage(); bool * blocked_dofs_val = blocked_dofs.storage(); UInt nb_degrees_of_freedom = x.getSize() * x.getNbComponent(); for (UInt n = 0; n < nb_degrees_of_freedom; ++n) { if(!(*blocked_dofs_val)) { *x_val = alpha * (*b_val / *A_val); } x_val++; A_val++; b_val++; blocked_dofs_val++; } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitPred() { AKANTU_DEBUG_IN(); if(increment_flag) { if(previous_displacement) increment->copy(*previous_displacement); else increment->copy(*displacement); } AKANTU_DEBUG_ASSERT(integrator,"itegrator should have been allocated: " << "have called initExplicit ? " << "or initImplicit ?"); integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *blocked_dofs); if(increment_flag) { Real * inc_val = increment->storage(); Real * dis_val = displacement->storage(); UInt nb_degree_of_freedom = displacement->getSize() * displacement->getNbComponent(); for (UInt n = 0; n < nb_degree_of_freedom; ++n) { *inc_val = *dis_val - *inc_val; inc_val++; dis_val++; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitCorr() { AKANTU_DEBUG_IN(); integrator->integrationSchemeCorrAccel(time_step, *displacement, *velocity, *acceleration, *blocked_dofs, *increment_acceleration); if(previous_displacement) previous_displacement->copy(*displacement); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStep() { AKANTU_DEBUG_IN(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeSolveStepEvent(method)); this->explicitPred(); this->updateResidual(); this->updateAcceleration(); this->explicitCorr(); EventManager::sendEvent(SolidMechanicsModelEvent::AfterSolveStepEvent(method)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Implicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * Initialize the solver and create the sparse matrices needed. * */ void SolidMechanicsModel::initSolver(__attribute__((unused)) SolverOptions & options) { #if !defined(AKANTU_USE_MUMPS) && !defined(AKANTU_USE_PETSC)// or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake AKANTU_DEBUG_ERROR("You should at least activate one solver."); #else UInt nb_global_nodes = mesh.getNbGlobalNodes(); delete jacobian_matrix; std::stringstream sstr; sstr << id << ":jacobian_matrix"; #ifdef AKANTU_USE_PETSC jacobian_matrix = new PETScMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id); #else jacobian_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id); #endif //AKANTU_USE PETSC jacobian_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension); if (!isExplicit()) { delete stiffness_matrix; std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix"; #ifdef AKANTU_USE_PETSC stiffness_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id); stiffness_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension); #else stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id); #endif //AKANTU_USE_PETSC } delete solver; std::stringstream sstr_solv; sstr_solv << id << ":solver"; #ifdef AKANTU_USE_PETSC solver = new SolverPETSc(*jacobian_matrix, sstr_solv.str()); #elif defined(AKANTU_USE_MUMPS) solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS if(solver) solver->initialize(options); #endif //AKANTU_HAS_SOLVER } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initJacobianMatrix() { #if defined(AKANTU_USE_MUMPS) && !defined(AKANTU_USE_PETSC) // @todo make it more flexible: this is an ugly patch to treat the case of non // fix profile of the K matrix delete jacobian_matrix; std::stringstream sstr_sti; sstr_sti << id << ":jacobian_matrix"; jacobian_matrix = new SparseMatrix(*stiffness_matrix, sstr_sti.str(), memory_id); std::stringstream sstr_solv; sstr_solv << id << ":solver"; delete solver; solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); if(solver) solver->initialize(_solver_no_options); #else AKANTU_DEBUG_ERROR("You need to activate the solver mumps."); #endif } /* -------------------------------------------------------------------------- */ /** * Initialize the implicit solver, either for dynamic or static cases, * * @param dynamic */ void SolidMechanicsModel::initImplicit(bool dynamic, SolverOptions & solver_options) { AKANTU_DEBUG_IN(); method = dynamic ? _implicit_dynamic : _static; if (!increment) setIncrementFlagOn(); initSolver(solver_options); if(method == _implicit_dynamic) { if(integrator) delete integrator; integrator = new TrapezoidalRule2(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initialAcceleration() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ma = f"); Solver * acc_solver = NULL; std::stringstream sstr; sstr << id << ":tmp_mass_matrix"; SparseMatrix * tmp_mass = new SparseMatrix(*mass_matrix, sstr.str(), memory_id); #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solver; sstr << id << ":solver_mass_matrix"; acc_solver = new SolverMumps(*mass_matrix, sstr_solver.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS acc_solver->initialize(); tmp_mass->applyBoundary(*blocked_dofs); acc_solver->setRHS(*residual); acc_solver->solve(*acceleration); delete acc_solver; delete tmp_mass; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); stiffness_matrix->clear(); // call compute stiffness matrix on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->assembleStiffnessMatrix(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseMatrix & SolidMechanicsModel::initVelocityDampingMatrix() { if(!velocity_damping_matrix) velocity_damping_matrix = new SparseMatrix(*jacobian_matrix, id + ":velocity_damping_matrix", memory_id); return *velocity_damping_matrix; } /* -------------------------------------------------------------------------- */ template<> bool SolidMechanicsModel::testConvergence<_scc_increment>(Real tolerance, Real & error){ AKANTU_DEBUG_IN(); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); error = 0; Real norm[2] = {0., 0.}; Real * increment_val = increment->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); Real * displacement_val = displacement->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); for (UInt d = 0; d < nb_degree_of_freedom; ++d) { if(!(*blocked_dofs_val) && is_local_node) { norm[0] += *increment_val * *increment_val; norm[1] += *displacement_val * *displacement_val; } blocked_dofs_val++; increment_val++; displacement_val++; } } StaticCommunicator::getStaticCommunicator().allReduce(norm, 2, _so_sum); norm[0] = sqrt(norm[0]); norm[1] = sqrt(norm[1]); AKANTU_DEBUG_ASSERT(!Math::isnan(norm[0]), "Something goes wrong in the solve phase"); if (norm[1] < Math::getTolerance()) { error = norm[0]; AKANTU_DEBUG_OUT(); // cout<<"Error 1: "< Math::getTolerance()) error = norm[0] / norm[1]; else error = norm[0]; //In case the total displacement is zero! // cout<<"Error 2: "< bool SolidMechanicsModel::testConvergence<_scc_residual>(Real tolerance, Real & norm) { AKANTU_DEBUG_IN(); UInt nb_nodes = residual->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); norm = 0; Real * residual_val = residual->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); if(is_local_node) { for (UInt d = 0; d < nb_degree_of_freedom; ++d) { if(!(*blocked_dofs_val)) { norm += *residual_val * *residual_val; } blocked_dofs_val++; residual_val++; } } else { blocked_dofs_val += spatial_dimension; residual_val += spatial_dimension; } } StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum); norm = sqrt(norm); AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (norm < tolerance); } /* -------------------------------------------------------------------------- */ template<> bool SolidMechanicsModel::testConvergence<_scc_residual_mass_wgh>(Real tolerance, Real & norm) { AKANTU_DEBUG_IN(); UInt nb_nodes = residual->getSize(); norm = 0; Real * residual_val = residual->storage(); Real * mass_val = this->mass->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); if(is_local_node) { for (UInt d = 0; d < spatial_dimension; ++d) { if(!(*blocked_dofs_val)) { norm += *residual_val * *residual_val/(*mass_val * *mass_val); } blocked_dofs_val++; residual_val++; mass_val++; } } else { blocked_dofs_val += spatial_dimension; residual_val += spatial_dimension; mass_val += spatial_dimension; } } StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum); norm = sqrt(norm); AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (norm < tolerance); } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance){ AKANTU_DEBUG_IN(); Real error=0; bool res = this->testConvergence<_scc_residual>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance, Real & error){ AKANTU_DEBUG_IN(); bool res = this->testConvergence<_scc_residual>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance){ AKANTU_DEBUG_IN(); Real error=0; bool res = this->testConvergence<_scc_increment>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance, Real & error){ AKANTU_DEBUG_IN(); bool res = this->testConvergence<_scc_increment>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitPred() { AKANTU_DEBUG_IN(); if(previous_displacement) previous_displacement->copy(*displacement); if(method == _implicit_dynamic) integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *blocked_dofs); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitCorr() { AKANTU_DEBUG_IN(); if(method == _implicit_dynamic) { integrator->integrationSchemeCorrDispl(time_step, *displacement, *velocity, *acceleration, *blocked_dofs, *increment); } else { UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; Real * incr_val = increment->storage(); Real * disp_val = displacement->storage(); bool * boun_val = blocked_dofs->storage(); for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, ++boun_val){ *incr_val *= (1. - *boun_val); *disp_val += *incr_val; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateIncrement() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(previous_displacement,"The previous displacement has to be initialized." << " Are you working with Finite or Ineslactic deformations?"); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; Real * incr_val = increment->storage(); Real * disp_val = displacement->storage(); Real * prev_disp_val = previous_displacement->storage(); for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, ++prev_disp_val) *incr_val = (*disp_val - *prev_disp_val); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updatePreviousDisplacement() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(previous_displacement,"The previous displacement has to be initialized." << " Are you working with Finite or Ineslactic deformations?"); previous_displacement->copy(*displacement); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeBoundaries() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized." << " Did you call initParallel?"); synch_registry->synchronize(_gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeResidual() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized." << " Did you call initPBC?"); synch_registry->synchronize(_gst_smm_res); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setIncrementFlagOn() { AKANTU_DEBUG_IN(); if(!increment) { UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_inc; sstr_inc << id << ":increment"; increment = &(alloc(sstr_inc.str(), nb_nodes, spatial_dimension, 0.)); } increment_flag = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real min_dt = getStableTimeStep(_not_ghost); /// reduction min over all processors StaticCommunicator::getStaticCommunicator().allReduce(&min_dt, 1, _so_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Material ** mat_val = &(materials.at(0)); Real min_dt = std::numeric_limits::max(); updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; elem.kind = _ek_regular; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type); for(; it != end; ++it) { elem.type = *it; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it); UInt nb_element = mesh.getNbElement(*it); Array::const_scalar_iterator mat_indexes = material_index(*it, ghost_type).begin(); Array::const_scalar_iterator mat_loc_num = material_local_numbering(*it, ghost_type).begin(); Array X(0, nb_nodes_per_element*spatial_dimension); FEEngine::extractNodalToElementField(mesh, *current_position, X, *it, _not_ghost); Array::matrix_iterator X_el = X.begin(spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++mat_indexes, ++mat_loc_num) { elem.element = *mat_loc_num; Real el_h = getFEEngine().getElementInradius(*X_el, *it); Real el_c = mat_val[*mat_indexes]->getCelerity(elem); Real el_dt = el_h / el_c; min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real energy = 0.; /// call update residual on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { energy += (*mat_it)->getPotentialEnergy(); } /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); if (!mass && !mass_matrix) AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); Real * vel_val = velocity->storage(); Real * mass_val = mass->storage(); for (UInt n = 0; n < nb_nodes; ++n) { Real mv2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !isPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < spatial_dimension; ++i) { if (count_node) mv2 += *vel_val * *vel_val * *mass_val; vel_val++; mass_val++; } ekin += mv2; } StaticCommunicator::getStaticCommunicator().allReduce(&ekin, 1, _so_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy(const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(type); Array vel_on_quad(nb_quadrature_points, spatial_dimension); Array filter_element(1, 1, index); getFEEngine().interpolateOnQuadraturePoints(*velocity, vel_on_quad, spatial_dimension, type, _not_ghost, filter_element); Array::vector_iterator vit = vel_on_quad.begin(spatial_dimension); Array::vector_iterator vend = vel_on_quad.end(spatial_dimension); Vector rho_v2(nb_quadrature_points); Real rho = materials[material_index(type)(index)]->getRho(); for (UInt q = 0; vit != vend; ++vit, ++q) { rho_v2(q) = rho * vit->dot(*vit); } AKANTU_DEBUG_OUT(); return .5*getFEEngine().integrate(rho_v2, type, index); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); Real * velo = velocity->storage(); Real * forc = force->storage(); Real * resi = residual->storage(); bool * boun = blocked_dofs->storage(); Real work = 0.; UInt nb_nodes = mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !isPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < spatial_dimension; ++i) { if (count_node) { if(*boun) work -= *resi * *velo * time_step; else work += *forc * *velo * time_step; } ++velo; ++forc; ++resi; ++boun; } } StaticCommunicator::getStaticCommunicator().allReduce(&work, 1, _so_sum); AKANTU_DEBUG_OUT(); return work; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(); } else if (energy_id == "external work"){ return getExternalWork(); } Real energy = 0.; std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { energy += (*mat_it)->getEnergy(energy_id); } /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id, const ElementType & type, UInt index){ AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(type, index); } std::vector::iterator mat_it; UInt mat_index = this->material_index(type, _not_ghost)(index); UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index); Real energy = this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); this->getFEEngine().initShapeFunctions(_not_ghost); this->getFEEngine().initShapeFunctions(_ghost); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = this->mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = this->mesh.lastType(spatial_dimension, gt, _ek_not_defined); for(; it != end; ++it) { UInt nb_element = this->mesh.getNbElement(*it, gt); if(!material_index.exists(*it, gt)) { this->material_index .alloc(nb_element, 1, *it, gt); this->material_local_numbering.alloc(nb_element, 1, *it, gt); } else { this->material_index (*it, gt).resize(nb_element); this->material_local_numbering(*it, gt).resize(nb_element); } } } Array::const_iterator it = element_list.begin(); Array::const_iterator end = element_list.end(); ElementTypeMapArray filter("new_element_filter", this->getID()); for (UInt el = 0; it != end; ++it, ++el) { const Element & elem = *it; if(!filter.exists(elem.type, elem.ghost_type)) filter.alloc(0, 1, elem.type, elem.ghost_type); filter(elem.type, elem.ghost_type).push_back(elem.element); } this->assignMaterialToElements(&filter); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsAdded(element_list, event); } if(method == _explicit_lumped_mass) this->assembleMassLumped(); if (method != _explicit_lumped_mass) { this->initSolver(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved(__attribute__((unused)) const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) { this->getFEEngine().initShapeFunctions(_not_ghost); this->getFEEngine().initShapeFunctions(_ghost); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsRemoved(element_list, new_numbering, event); } - if (method != _explicit_lumped_mass) { - this->initSolver(); - } - } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Array & nodes_list, __attribute__((unused)) const NewNodesEvent & event) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); if(displacement) displacement->resize(nb_nodes); if(mass ) mass ->resize(nb_nodes); if(velocity ) velocity ->resize(nb_nodes); if(acceleration) acceleration->resize(nb_nodes); if(force ) force ->resize(nb_nodes); if(residual ) residual ->resize(nb_nodes); if(blocked_dofs) blocked_dofs->resize(nb_nodes); if(previous_displacement) previous_displacement->resize(nb_nodes); if(increment_acceleration) increment_acceleration->resize(nb_nodes); if(increment) increment->resize(nb_nodes); if(current_position) current_position->resize(nb_nodes); delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onNodesAdded(nodes_list, event); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesRemoved(__attribute__((unused)) const Array & element_list, const Array & new_numbering, __attribute__((unused)) const RemovedNodesEvent & event) { if(displacement) mesh.removeNodesFromArray(*displacement, new_numbering); if(mass ) mesh.removeNodesFromArray(*mass , new_numbering); if(velocity ) mesh.removeNodesFromArray(*velocity , new_numbering); if(acceleration) mesh.removeNodesFromArray(*acceleration, new_numbering); if(force ) mesh.removeNodesFromArray(*force , new_numbering); if(residual ) mesh.removeNodesFromArray(*residual , new_numbering); if(blocked_dofs) mesh.removeNodesFromArray(*blocked_dofs, new_numbering); if(increment_acceleration) mesh.removeNodesFromArray(*increment_acceleration, new_numbering); if(increment) mesh.removeNodesFromArray(*increment , new_numbering); delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); + + if (method != _explicit_lumped_mass) { + this->initSolver(); + } + } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::isInternal(const std::string & field_name, const ElementKind & element_kind){ bool is_internal = false; /// check if at least one material contains field_id as an internal for (UInt m = 0; m < materials.size() && !is_internal; ++m) { is_internal |= materials[m]->isInternal(field_name, element_kind); } return is_internal; } /* -------------------------------------------------------------------------- */ ElementTypeMap SolidMechanicsModel::getInternalDataPerElem(const std::string & field_name, const ElementKind & element_kind){ if (!(this->isInternal(field_name,element_kind))) AKANTU_EXCEPTION("unknown internal " << field_name); for (UInt m = 0; m < materials.size() ; ++m) { if (materials[m]->isInternal(field_name, element_kind)) return materials[m]->getInternalDataPerElem(field_name, element_kind); } return ElementTypeMap(); } /* -------------------------------------------------------------------------- */ ElementTypeMapArray & SolidMechanicsModel::flattenInternal(const std::string & field_name, const ElementKind & kind, const GhostType ghost_type){ std::pair key(field_name,kind); if (this->registered_internals.count(key) == 0){ this->registered_internals[key] = new ElementTypeMapArray(field_name, this->id); } ElementTypeMapArray * internal_flat = this->registered_internals[key]; + + typedef ElementTypeMapArray::type_iterator iterator; + iterator tit = internal_flat->firstType(spatial_dimension, + ghost_type, + kind); + iterator end = internal_flat->lastType(spatial_dimension, + ghost_type, + kind); + + for (; tit != end; ++tit) { + ElementType type = *tit; + (*internal_flat)(type,ghost_type).clear(); + } + for (UInt m = 0; m < materials.size(); ++m) { if (materials[m]->isInternal(field_name, kind)) materials[m]->flattenInternal(field_name, *internal_flat, ghost_type, kind); } return *internal_flat; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::flattenAllRegisteredInternals(const ElementKind & kind){ + std::map, ElementTypeMapArray *>::iterator it = this->registered_internals.begin(); std::map, ElementTypeMapArray *>::iterator end = this->registered_internals.end(); while (it != end){ if (kind == it->first.second) this->flattenInternal(it->first.first, kind); ++it; } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onDump(){ this->flattenAllRegisteredInternals(_ek_regular); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER dumper::Field * SolidMechanicsModel ::createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) { dumper::Field * field = NULL; if(field_name == "partitions") field = mesh.createElementalField(mesh.getConnectivities(),group_name,spatial_dimension,kind); else if(field_name == "material_index") field = mesh.createElementalField(material_index,group_name,spatial_dimension,kind); else { // this copy of field_name is used to compute derivated data such as // strain and von mises stress that are based on grad_u and stress std::string field_name_copy(field_name); if (field_name == "strain" || field_name == "Green strain" || field_name == "principal strain" || field_name == "principal Green strain") field_name_copy = "grad_u"; else if (field_name == "Von Mises stress") field_name_copy = "stress"; bool is_internal = this->isInternal(field_name_copy,kind); if (is_internal) { ElementTypeMap nb_data_per_elem = this->getInternalDataPerElem(field_name_copy, kind); ElementTypeMapArray & internal_flat = this->flattenInternal(field_name_copy,kind); field = mesh.createElementalField(internal_flat, group_name, spatial_dimension,kind,nb_data_per_elem); if (field_name == "strain"){ dumper::ComputeStrain * foo = new dumper::ComputeStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } else if (field_name == "Von Mises stress") { dumper::ComputeVonMisesStress * foo = new dumper::ComputeVonMisesStress(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } else if (field_name == "Green strain") { dumper::ComputeStrain * foo = new dumper::ComputeStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } else if (field_name == "principal strain") { dumper::ComputePrincipalStrain * foo = new dumper::ComputePrincipalStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } else if (field_name == "principal Green strain") { dumper::ComputePrincipalStrain * foo = new dumper::ComputePrincipalStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } //treat the paddings if (padding_flag){ if (field_name == "stress"){ if (spatial_dimension == 2) { dumper::StressPadder<2> * foo = new dumper::StressPadder<2>(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } } else if (field_name == "strain" || field_name == "Green strain"){ if (spatial_dimension == 2) { dumper::StrainPadder<2> * foo = new dumper::StrainPadder<2>(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } } } // homogenize the field dumper::ComputeFunctorInterface * foo = dumper::HomogenizerProxy::createHomogenizer(*field); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } } return field; } /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { std::map* > real_nodal_fields; real_nodal_fields["displacement" ] = displacement; real_nodal_fields["mass" ] = mass; real_nodal_fields["velocity" ] = velocity; real_nodal_fields["acceleration" ] = acceleration; real_nodal_fields["force" ] = force; real_nodal_fields["residual" ] = residual; real_nodal_fields["increment" ] = increment; dumper::Field * field = NULL; if (padding_flag) field = mesh.createNodalField(real_nodal_fields[field_name],group_name, 3); else field = mesh.createNodalField(real_nodal_fields[field_name],group_name); return field; } /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel::createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) { std::map* > uint_nodal_fields; uint_nodal_fields["blocked_dofs" ] = blocked_dofs; dumper::Field * field = NULL; field = mesh.createNodalField(uint_nodal_fields[field_name],group_name); return field; } /* -------------------------------------------------------------------------- */ #else /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel ::createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind){ return NULL; } /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { return NULL; } /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel::createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) { return NULL; } #endif /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); synch_registry->synchronize(_gst_for_dump); mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name, UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); synch_registry->synchronize(_gst_for_dump); mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name, Real time, UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); synch_registry->synchronize(_gst_for_dump); mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump() { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(Real time, UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(time, step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeCauchyStresses() { AKANTU_DEBUG_IN(); // call compute stiffness matrix on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; if(mat.isFiniteDeformation()) mat.computeAllCauchyStresses(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::saveStressAndStrainBeforeDamage() { EventManager::sendEvent(SolidMechanicsModelEvent::BeginningOfDamageIterationEvent()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateEnergiesAfterDamage() { EventManager::sendEvent(SolidMechanicsModelEvent::AfterDamageEvent()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Solid Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); mass ->printself(stream, indent + 2); velocity ->printself(stream, indent + 2); acceleration->printself(stream, indent + 2); force ->printself(stream, indent + 2); residual ->printself(stream, indent + 2); blocked_dofs->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + material information [" << std::endl; material_index.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + materials [" << std::endl; std::vector::const_iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { const Material & mat = *(*mat_it); mat.printself(stream, indent + 1); } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__