diff --git a/src/model/solid_mechanics/material.cc b/src/model/solid_mechanics/material.cc index 5ec509076..270238636 100644 --- a/src/model/solid_mechanics/material.cc +++ b/src/model/solid_mechanics/material.cc @@ -1,1376 +1,1378 @@ /** * @file material.cc * * @author Aurelia Isabel Cuba Ramos * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Nov 24 2015 * * @brief Implementation of the common part of the material class * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 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 . * */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, const ID & id) : Memory(id, model.getMemoryID()), Parsable(ParserType::_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), eigengradu("eigen_grad_u", *this), gradu("grad_u", *this), green_strain("green_strain", *this), piola_kirchhoff_2("piola_kirchhoff_2", *this), 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 element_filter.initialize(model.getMesh(), _spatial_dimension = spatial_dimension); // 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(ParserType::_material, id), is_init(false), fem(fe_engine), 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), eigengradu("eigen_grad_u", *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("piola_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(); element_filter.initialize(mesh, _spatial_dimension = spatial_dimension); // mesh.initElementTypeMapArray(element_filter, 1, spatial_dimension, false, // _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::~Material() = default; /* -------------------------------------------------------------------------- */ void Material::initialize() { registerParam("rho", rho, Real(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 eigengradu.initialize(spatial_dimension * spatial_dimension); gradu.initialize(spatial_dimension * spatial_dimension); stress.initialize(spatial_dimension * spatial_dimension); potential_energy.initialize(1); 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 (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->resize(); for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->resize(); for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->resize(); is_init = true; + updateInternalParameters(); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::savePreviousState() { AKANTU_DEBUG_IN(); for (auto 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::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); if (!finite_deformation) { auto & internal_force = const_cast &>(model.getInternalForce()); // Mesh & mesh = fem.getMesh(); for (auto && type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0) continue; const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(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"); fem.computeBtD(stress(type, ghost_type), *sigma_dphi_dx, type, ghost_type, elem_filter); /** * 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, type, ghost_type, elem_filter); delete sigma_dphi_dx; /// assemble model.getDOFManager().assembleElementalArrayLocalArray( *int_sigma_dphi_dx, internal_force, type, ghost_type, -1, elem_filter); delete int_sigma_dphi_dx; } } else { switch (spatial_dimension) { case 1: this->assembleInternalForces<1>(ghost_type); break; case 2: this->assembleInternalForces<2>(ghost_type); break; case 3: this->assembleInternalForces<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(); for (const auto & type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.size() == 0) continue; Array & gradu_vect = gradu(type, ghost_type); /// compute @f$\nabla u@f$ fem.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vect, spatial_dimension, type, ghost_type, elem_filter); gradu_vect -= eigengradu(type, ghost_type); /// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$ computeStress(type, 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."); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { switch (spatial_dimension) { case 1: this->computeCauchyStress<1>(type, ghost_type); break; case 2: this->computeCauchyStress<2>(type, ghost_type); break; case 3: this->computeCauchyStress<3>(type, 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(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); /// compute @f$\nabla u@f$ fem.gradientOnIntegrationPoints(displacement, gradu_vect, spatial_dimension, type, ghost_type, elem_filter); setToSteadyState(type, 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(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { if (finite_deformation) { switch (spatial_dimension) { case 1: { assembleStiffnessMatrixNL<1>(type, ghost_type); assembleStiffnessMatrixL2<1>(type, ghost_type); break; } case 2: { assembleStiffnessMatrixNL<2>(type, ghost_type); assembleStiffnessMatrixL2<2>(type, ghost_type); break; } case 3: { assembleStiffnessMatrixNL<3>(type, ghost_type); assembleStiffnessMatrixL2<3>(type, ghost_type); break; } } } else { switch (spatial_dimension) { case 1: { assembleStiffnessMatrix<1>(type, ghost_type); break; } case 2: { assembleStiffnessMatrix<2>(type, ghost_type); break; } case 3: { assembleStiffnessMatrix<3>(type, 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.size() == 0) { AKANTU_DEBUG_OUT(); return; } // const Array & shapes_derivatives = // fem.getShapesDerivatives(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); fem.gradientOnIntegrationPoints(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); /// 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"); fem.computeBtDB(*tangent_stiffness_matrix, *bt_d_b, 4, type, ghost_type, elem_filter); delete tangent_stiffness_matrix; /// 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; model.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixNL(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 = delta_gradu(type, ghost_type); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); Array * shapes_derivatives_filtered = new Array( nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); fem.filterElementalData(fem.getMesh(), shapes_derivatives, *shapes_derivatives_filtered, type, ghost_type, elem_filter); /// 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); auto shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin( spatial_dimension, nb_nodes_per_element); auto Bt_S_B_it = bt_s_b->begin(bt_s_b_size, bt_s_b_size); auto Bt_S_B_end = bt_s_b->end(bt_s_b_size, bt_s_b_size); auto 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) { auto & Bt_S_B = *Bt_S_B_it; const auto & Piola_kirchhoff_matrix = *piola_it; setCauchyStressMatrix(Piola_kirchhoff_matrix, S); VoigtHelper::transferBMatrixToBNL(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bt_S.template mul(B, S); Bt_S_B.template mul(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; model.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixL2(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.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); fem.gradientOnIntegrationPoints(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"); fem.filterElementalData(fem.getMesh(), shapes_derivatives, *shapes_derivatives_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 B2(tangent_size, dim * nb_nodes_per_element); Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); auto shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin( spatial_dimension, nb_nodes_per_element); auto Bt_D_B_it = bt_d_b->begin(bt_d_b_size, bt_d_b_size); auto grad_u_it = gradu_vect.begin(dim, dim); auto D_it = tangent_stiffness_matrix->begin(tangent_size, tangent_size); auto 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) { const auto & grad_u = *grad_u_it; const auto & D = *D_it; auto & 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.template mul(B, D); Bt_D_B.template mul(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; model.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & internal_force = model.getInternalForce(); Mesh & mesh = fem.getMesh(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.size() == 0) continue; UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); Array * shapesd_filtered = new Array( nb_element, size_of_shapes_derivatives, "filtered shapesd"); fem.filterElementalData(mesh, shapes_derivatives, *shapesd_filtered, type, 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; auto * bt_s = new Array(nb_element * nb_quadrature_points, bt_s_size, "B^t*S"); auto grad_u_it = this->gradu(type, ghost_type).begin(dim, dim); auto grad_u_end = this->gradu(type, ghost_type).end(dim, dim); auto stress_it = this->piola_kirchhoff_2(type, 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) { auto & grad_u = *grad_u_it; auto & r_it = *bt_s_it; auto & 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.template mul(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, type, ghost_type, elem_filter); delete bt_s; model.getDOFManager().assembleElementalArrayLocalArray( *r_e, internal_force, type, ghost_type, -1., elem_filter); delete r_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergyByElements() { AKANTU_DEBUG_IN(); for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { computePotentialEnergy(type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergy(ElementType, GhostType) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real epot = 0.; computePotentialEnergyByElements(); /// integrate the potential energy for each type of elements for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { epot += fem.integrate(potential_energy(type, _not_ghost), type, _not_ghost, element_filter(type, _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.getNbIntegrationPoints(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(const std::string & type) { AKANTU_DEBUG_IN(); if (type == "potential") return getPotentialEnergy(); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(const 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::initElementalFieldInterpolation( const ElementTypeMapArray & interpolation_points_coordinates) { AKANTU_DEBUG_IN(); this->fem.initElementalFieldInterpolationFromIntegrationPoints( interpolation_points_coordinates, this->interpolation_points_matrices, this->interpolation_inverse_coordinates, &(this->element_filter)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::interpolateStress(ElementTypeMapArray & result, const GhostType ghost_type) { this->fem.interpolateElementalFieldFromIntegrationPoints( this->stress, this->interpolation_points_matrices, this->interpolation_inverse_coordinates, result, ghost_type, &(this->element_filter)); } /* -------------------------------------------------------------------------- */ void Material::interpolateStressOnFacets( ElementTypeMapArray & result, ElementTypeMapArray & by_elem_result, const GhostType ghost_type) { interpolateStress(by_elem_result, ghost_type); UInt stress_size = this->stress.getNbComponent(); const Mesh & mesh = this->model.getMesh(); const Mesh & mesh_facets = mesh.getMeshFacets(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_fil = element_filter(type, ghost_type); Array & by_elem_res = by_elem_result(type, ghost_type); UInt nb_element = elem_fil.size(); UInt nb_element_full = this->model.getMesh().getNbElement(type, ghost_type); UInt nb_interpolation_points_per_elem = by_elem_res.size() / nb_element_full; 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> * element_to_facet = nullptr; GhostType current_ghost_type = _casper; Array * result_vec = nullptr; Array::const_matrix_iterator result_it = by_elem_res.begin_reinterpret( stress_size, nb_interpolation_points_per_elem, nb_element_full); for (UInt el = 0; el < nb_element; ++el) { 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 * stress_size, stress_size); const Matrix & result_tmp(result_it[global_el]); result_local = result_tmp(f * nb_quad_per_facet + q); } } } } } /* -------------------------------------------------------------------------- */ template const Array & Material::getArray(const ID & /*vect_id*/, const ElementType & /*type*/, const GhostType & /*ghost_type*/) const { AKANTU_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template Array & Material::getArray(const ID & /*vect_id*/, const ElementType & /*type*/, const GhostType & /*ghost_type*/) { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ 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 << "(" << getID() << ") does not contain a vector " << vect_id << " (" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template <> 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 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 << "(" << getID() << ") does not contain a vector " << vect_id << " (" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template <> 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(__attribute__((unused)) const ID & int_id) const { AKANTU_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template InternalField & Material::getInternal(__attribute__((unused)) const ID & int_id) { AKANTU_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template <> const InternalField & Material::getInternal(const ID & int_id) const { auto 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) { auto 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 <> const InternalField & Material::getInternal(const ID & int_id) const { auto 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; } /* -------------------------------------------------------------------------- */ template <> InternalField & Material::getInternal(const ID & int_id) { auto 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(), getMemoryID()); 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.size(), 1, type, ghost_type); Array & mat_renumbering = material_local_new_numbering(type, ghost_type); UInt nb_element = elem_filter.size(); 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.size()); elem_filter.copy(elem_filter_tmp); } } for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::resizeInternals() { AKANTU_DEBUG_IN(); for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->resize(); for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->resize(); for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->resize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::onElementsAdded(const Array &, const NewElementsEvent &) { 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(), getMemoryID()); auto el_begin = element_list.begin(); auto el_end = element_list.end(); for (auto && gt : ghost_types) { for (auto && type : new_numbering.elementTypes(_all_dimensions, gt, _ek_not_defined)) { if (not element_filter.exists(type, gt) || element_filter(type, gt).size() == 0) continue; auto & elem_filter = element_filter(type, gt); auto & mat_indexes = this->model.getMaterialByElement(type, gt); auto & mat_loc_num = this->model.getMaterialLocalNumbering(type, gt); auto 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.size(), 1, type, gt); auto & mat_renumbering = material_local_new_numbering(type, gt); const auto & renumbering = new_numbering(type, gt); Array elem_filter_tmp; UInt ni = 0; Element el{type, 0, gt}; for (UInt i = 0; i < elem_filter.size(); ++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.size()); elem_filter.copy(elem_filter_tmp); } } for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); } /* -------------------------------------------------------------------------- */ void Material::beforeSolveStep() { this->savePreviousState(); } /* -------------------------------------------------------------------------- */ void Material::afterSolveStep() { for (auto & type : element_filter.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { this->updateEnergies(type, _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) { 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; } /* -------------------------------------------------------------------------- */ /// extrapolate internal values void Material::extrapolateInternal(const ID & id, const Element & element, __attribute__((unused)) const Matrix & point, Matrix & extrapolated) { if (this->isInternal(id, element.kind())) { UInt nb_element = this->element_filter(element.type, element.ghost_type).size(); const ID name = this->getID() + ":" + id; UInt nb_quads = this->internal_vectors_real[name]->getFEEngine().getNbIntegrationPoints( 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!!!! AKANTU_DEBUG_WARNING("This is a fix, values are not truly extrapolated"); 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; break; } } for (UInt i = 0; i < extrapolated.size(); ++i) { extrapolated(i) = values(index); } } else { Matrix default_values(extrapolated.rows(), extrapolated.cols(), 0.); extrapolated = default_values; } } /* -------------------------------------------------------------------------- */ void Material::applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const GhostType ghost_type) { for (auto && type : element_filter.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { if (!element_filter(type, ghost_type).size()) continue; auto eigen_it = this->eigengradu(type, ghost_type) .begin(spatial_dimension, spatial_dimension); auto eigen_end = this->eigengradu(type, ghost_type) .end(spatial_dimension, spatial_dimension); for (; eigen_it != eigen_end; ++eigen_it) { auto & current_eigengradu = *eigen_it; current_eigengradu = prescribed_eigen_grad_u; } } } } // namespace akantu diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh index c6b70452a..3a73163a2 100644 --- a/src/model/solid_mechanics/material.hh +++ b/src/model/solid_mechanics/material.hh @@ -1,684 +1,683 @@ /** * @file material.hh * * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Wed Nov 25 2015 * * @brief Mother class for all materials * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_factory.hh" #include "aka_voigthelper.hh" #include "data_accessor.hh" #include "integration_point.hh" #include "parsable.hh" #include "parser.hh" /* -------------------------------------------------------------------------- */ #include "internal_field.hh" #include "random_internal_field.hh" /* -------------------------------------------------------------------------- */ #include "mesh_events.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_HH__ #define __AKANTU_MATERIAL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class SolidMechanicsModel; } // namespace akantu namespace 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: #if __cplusplus > 199711L Material(const Material & mat) = delete; Material & operator=(const Material & mat) = delete; #endif /// 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 ~Material() override; 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_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_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_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(__attribute__((unused)) const Element & element) const { AKANTU_TO_IMPLEMENT(); } /// compute the s-wave speed in the material virtual Real getShearWaveSpeed(__attribute__((unused)) const Element & element) const { AKANTU_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_TO_IMPLEMENT(); } template void unregisterInternal(__attribute__((unused)) InternalField & vect) { AKANTU_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 assembleInternalForces(GhostType ghost_type); /// 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 // 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); inline UInt addElement(const Element & element); /// 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 void printself(std::ostream & stream, int indent = 0) const override; /** * 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, ElementTypeMapArray & by_elem_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: /* ------------------------------------------------------------------------ */ inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const; /// compute the potential energy by element void computePotentialEnergyByElements(); /// resize the intenals arrays virtual void resizeInternals(); /* ------------------------------------------------------------------------ */ /* Finite deformation functions */ /* This functions area implementing what is described in the paper of Bathe */ /* et al, in IJNME, Finite Element Formulations For Large Deformation */ /* Dynamic Analysis, Vol 9, 353-386, 1975 */ /* ------------------------------------------------------------------------ */ protected: /// assemble the residual template void assembleInternalForces(GhostType ghost_type); /// Computation of Cauchy stress tensor in the case of finite deformation from /// the 2nd Piola-Kirchhoff for a given element type template void computeCauchyStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// Computation the Cauchy stress the 2nd Piola-Kirchhoff and the deformation /// gradient 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); /// 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 & sigma); /// write the stress tensor in the Voigt notation. template inline void setCauchyStressArray(const Matrix & S_t, Matrix & sigma_voight); /* ------------------------------------------------------------------------ */ /* Conversion functions */ /* ------------------------------------------------------------------------ */ public: 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 IntegrationPoint convertToLocalPoint(const IntegrationPoint & global_point) const; /// converts local quadrature point to global quadrature point inline IntegrationPoint convertToGlobalPoint(const IntegrationPoint & local_point) const; /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; 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: /* ------------------------------------------------------------------------ */ void onNodesAdded(const Array &, const NewNodesEvent &) override{}; void onNodesRemoved(const Array &, const Array &, const RemovedNodesEvent &) override{}; void onElementsAdded(const Array & element_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array &, const Array &, const ElementTypeMapArray &, const ChangedElementsEvent &) override{}; /* ------------------------------------------------------------------------ */ /* SolidMechanicsModelEventHandler inherited members */ /* ------------------------------------------------------------------------ */ public: virtual void beforeSolveStep(); virtual void afterSolveStep(); void onDamageIteration() override; void onDamageUpdate() override; void onDump() override; /* ------------------------------------------------------------------------ */ /* 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(const std::string & energy_id); /// return the energy (identified by id) for the provided element virtual Real getEnergy(const 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; } 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); template const InternalField & getInternal(const ID & id) const; template InternalField & getInternal(const ID & id); template inline bool isInternal(const ID & id, const ElementKind & element_kind) const; template ElementTypeMap getInternalDataPerElem(const ID & id, const ElementKind & element_kind) const; bool isFiniteDeformation() const { return finite_deformation; } bool isInelasticDeformation() const { return inelastic_deformation; } template inline void setParam(const ID & param, T value); - template inline void setParamNoUpdate(const ID & param, T value); inline const Parameter & getParam(const ID & param) const; template void flattenInternal(const std::string & field_id, ElementTypeMapArray & internal_flat, const GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined) const; /// apply a constant eigengrad_u everywhere in the material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const GhostType = _not_ghost); /// specify if the matrix need to be recomputed for this material virtual bool hasStiffnessMatrixChanged() { return true; } protected: bool isInit() const { return is_init; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ 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; /// eigengrad_u arrays ordered by element types InternalField eigengradu; /// 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; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Material & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "material_inline_impl.cc" #include "internal_field_tmpl.hh" #include "random_internal_field_tmpl.hh" /* -------------------------------------------------------------------------- */ /* Auto loop */ /* -------------------------------------------------------------------------- */ /// This can be used to automatically write the loop on quadrature points in /// functions such as computeStress. This macro in addition to write the loop /// provides two tensors (matrices) sigma and grad_u #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).size()); \ \ Array::iterator> 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).size()); \ 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 } /// This can be used to automatically write the loop on quadrature points in /// functions such as computeTangentModuli. This macro in addition to write the /// loop provides two tensors (matrices) sigma_tensor, grad_u, and a matrix /// where the elemental tangent moduli should be stored in Voigt Notation #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).size()); \ \ 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 } /* -------------------------------------------------------------------------- */ namespace akantu { using MaterialFactory = Factory; } // namespace akantu #define INSTANTIATE_MATERIAL_ONLY(mat_name) \ template class mat_name<1>; \ template class mat_name<2>; \ template class mat_name<3> #define MATERIAL_DEFAULT_PER_DIM_ALLOCATOR(id, mat_name) \ [](UInt dim, const ID &, SolidMechanicsModel & model, \ const ID & id) -> std::unique_ptr { \ switch (dim) { \ case 1: \ return std::make_unique>(model, id); \ case 2: \ return std::make_unique>(model, id); \ case 3: \ return std::make_unique>(model, id); \ default: \ AKANTU_EXCEPTION("The dimension " \ << dim << "is not a valid dimension for the material " \ << #id); \ } \ } #define INSTANTIATE_MATERIAL(id, mat_name) \ INSTANTIATE_MATERIAL_ONLY(mat_name); \ static bool material_is_alocated_##id [[gnu::unused]] = \ MaterialFactory::getInstance().registerAllocator( \ #id, MATERIAL_DEFAULT_PER_DIM_ALLOCATOR(id, mat_name)) #endif /* __AKANTU_MATERIAL_HH__ */ diff --git a/src/model/solid_mechanics/material_inline_impl.cc b/src/model/solid_mechanics/material_inline_impl.cc index 303d71707..aaedb2e96 100644 --- a/src/model/solid_mechanics/material_inline_impl.cc +++ b/src/model/solid_mechanics/material_inline_impl.cc @@ -1,475 +1,464 @@ /** * @file material_inline_impl.cc * * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue Jul 27 2010 * @date last modification: Wed Nov 25 2015 * * @brief Implementation of the inline functions of the class material * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 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 . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_INLINE_IMPL_CC__ #define __AKANTU_MATERIAL_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ inline UInt Material::addElement(const ElementType & type, UInt element, const GhostType & ghost_type) { Array & el_filter = this->element_filter(type, ghost_type); el_filter.push_back(element); return el_filter.size() - 1; } /* -------------------------------------------------------------------------- */ inline UInt Material::addElement(const Element & element) { return this->addElement(element.type, element.element, element.ghost_type); } /* -------------------------------------------------------------------------- */ inline UInt Material::getTangentStiffnessVoigtSize(UInt dim) const { return (dim * (dim - 1) / 2 + dim); } /* -------------------------------------------------------------------------- */ inline UInt Material::getCauchyStressMatrixSize(UInt dim) const { return (dim * dim); } /* -------------------------------------------------------------------------- */ template inline void Material::gradUToF(const Matrix & grad_u, Matrix & F) { AKANTU_DEBUG_ASSERT(F.size() >= grad_u.size() && grad_u.size() == dim * dim, "The dimension of the tensor F should be greater or " "equal to the dimension of the tensor grad_u."); F.eye(); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) F(i, j) += grad_u(i, j); } /* -------------------------------------------------------------------------- */ template inline void Material::computeCauchyStressOnQuad(const Matrix & F, const Matrix & piola, Matrix & sigma, const Real & C33) const { Real J = F.det() * sqrt(C33); Matrix F_S(dim, dim); F_S.mul(F, piola); Real constant = J ? 1. / J : 0; sigma.mul(F_S, F, constant); } /* -------------------------------------------------------------------------- */ inline void Material::rightCauchy(const Matrix & F, Matrix & C) { C.mul(F, F); } /* -------------------------------------------------------------------------- */ inline void Material::leftCauchy(const Matrix & F, Matrix & B) { B.mul(F, F); } /* -------------------------------------------------------------------------- */ template inline void Material::gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon) { for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) epsilon(i, j) = 0.5 * (grad_u(i, j) + grad_u(j, i)); } /* -------------------------------------------------------------------------- */ template inline void Material::gradUToGreenStrain(const Matrix & grad_u, Matrix & epsilon) { epsilon.mul(grad_u, grad_u, .5); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) epsilon(i, j) += 0.5 * (grad_u(i, j) + grad_u(j, i)); } /* -------------------------------------------------------------------------- */ inline Real Material::stressToVonMises(const Matrix & stress) { // compute deviatoric stress UInt dim = stress.cols(); Matrix deviatoric_stress = Matrix::eye(dim, -1. * stress.trace() / 3.); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) deviatoric_stress(i, j) += stress(i, j); // return Von Mises stress return std::sqrt(3. * deviatoric_stress.doubleDot(deviatoric_stress) / 2.); } /* ---------------------------------------------------------------------------*/ template inline void Material::setCauchyStressArray(const Matrix & S_t, Matrix & sigma_voight) { AKANTU_DEBUG_IN(); sigma_voight.clear(); // see Finite element formulations for large deformation dynamic analysis, // Bathe et al. IJNME vol 9, 1975, page 364 ^t\tau /* * 1d: [ s11 ]' * 2d: [ s11 s22 s12 ]' * 3d: [ s11 s22 s33 s23 s13 s12 ] */ for (UInt i = 0; i < dim; ++i) // diagonal terms sigma_voight(i, 0) = S_t(i, i); for (UInt i = 1; i < dim; ++i) // term s12 in 2D and terms s23 s13 in 3D sigma_voight(dim + i - 1, 0) = S_t(dim - i - 1, dim - 1); for (UInt i = 2; i < dim; ++i) // term s13 in 3D sigma_voight(dim + i, 0) = S_t(0, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void Material::setCauchyStressMatrix(const Matrix & S_t, Matrix & sigma) { AKANTU_DEBUG_IN(); sigma.clear(); /// see Finite ekement formulations for large deformation dynamic analysis, /// Bathe et al. IJNME vol 9, 1975, page 364 ^t\tau for (UInt i = 0; i < dim; ++i) { for (UInt m = 0; m < dim; ++m) { for (UInt n = 0; n < dim; ++n) { sigma(i * dim + m, i * dim + n) = S_t(m, n); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline Element Material::convertToLocalElement(const Element & global_element) const { UInt ge = global_element.element; #ifndef AKANTU_NDEBUG UInt model_mat_index = this->model.getMaterialByElement( global_element.type, global_element.ghost_type)(ge); UInt mat_index = this->model.getMaterialIndex(this->name); AKANTU_DEBUG_ASSERT(model_mat_index == mat_index, "Conversion of a global element in a local element for " "the wrong material " << this->name << std::endl); #endif UInt le = this->model.getMaterialLocalNumbering( global_element.type, global_element.ghost_type)(ge); Element tmp_quad{global_element.type, le, global_element.ghost_type}; return tmp_quad; } /* -------------------------------------------------------------------------- */ inline Element Material::convertToGlobalElement(const Element & local_element) const { UInt le = local_element.element; UInt ge = this->element_filter(local_element.type, local_element.ghost_type)(le); Element tmp_quad{local_element.type, ge, local_element.ghost_type}; return tmp_quad; } /* -------------------------------------------------------------------------- */ inline IntegrationPoint Material::convertToLocalPoint(const IntegrationPoint & global_point) const { const FEEngine & fem = this->model.getFEEngine(); UInt nb_quad = fem.getNbIntegrationPoints(global_point.type); Element el = this->convertToLocalElement(static_cast(global_point)); IntegrationPoint tmp_quad(el, global_point.num_point, nb_quad); return tmp_quad; } /* -------------------------------------------------------------------------- */ inline IntegrationPoint Material::convertToGlobalPoint(const IntegrationPoint & local_point) const { const FEEngine & fem = this->model.getFEEngine(); UInt nb_quad = fem.getNbIntegrationPoints(local_point.type); Element el = this->convertToGlobalElement(static_cast(local_point)); IntegrationPoint tmp_quad(el, local_point.num_point, nb_quad); return tmp_quad; } /* -------------------------------------------------------------------------- */ inline UInt Material::getNbData(const Array & elements, const SynchronizationTag & tag) const { if (tag == _gst_smm_stress) { return (this->isFiniteDeformation() ? 3 : 1) * spatial_dimension * spatial_dimension * sizeof(Real) * this->getModel().getNbIntegrationPoints(elements); } return 0; } /* -------------------------------------------------------------------------- */ inline void Material::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { if (tag == _gst_smm_stress) { if (this->isFiniteDeformation()) { packElementDataHelper(piola_kirchhoff_2, buffer, elements); packElementDataHelper(gradu, buffer, elements); } packElementDataHelper(stress, buffer, elements); } } /* -------------------------------------------------------------------------- */ inline void Material::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { if (tag == _gst_smm_stress) { if (this->isFiniteDeformation()) { unpackElementDataHelper(piola_kirchhoff_2, buffer, elements); unpackElementDataHelper(gradu, buffer, elements); } unpackElementDataHelper(stress, buffer, elements); } } /* -------------------------------------------------------------------------- */ inline const Parameter & Material::getParam(const ID & param) const { try { return get(param); } catch (...) { AKANTU_EXCEPTION("No parameter " << param << " in the material " << getID()); } } /* -------------------------------------------------------------------------- */ template inline void Material::setParam(const ID & param, T value) { try { set(param, value); } catch (...) { AKANTU_EXCEPTION("No parameter " << param << " in the material " << getID()); } updateInternalParameters(); } -/* -------------------------------------------------------------------------- */ -template -inline void Material::setParamNoUpdate(const ID & param, T value) { - try { - set(param, value); - } catch (...) { - AKANTU_EXCEPTION("No parameter " << param << " in the material " - << getID()); - } -} - /* -------------------------------------------------------------------------- */ template inline void Material::packElementDataHelper( const ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id) const { DataAccessor::packElementalDataHelper(data_to_pack, buffer, elements, true, model.getFEEngine(fem_id)); } /* -------------------------------------------------------------------------- */ template inline void Material::unpackElementDataHelper( ElementTypeMapArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id) { DataAccessor::unpackElementalDataHelper(data_to_unpack, buffer, elements, true, model.getFEEngine(fem_id)); } /* -------------------------------------------------------------------------- */ template <> inline void Material::registerInternal(InternalField & vect) { internal_vectors_real[vect.getID()] = &vect; } template <> inline void Material::registerInternal(InternalField & vect) { internal_vectors_uint[vect.getID()] = &vect; } template <> inline void Material::registerInternal(InternalField & vect) { internal_vectors_bool[vect.getID()] = &vect; } /* -------------------------------------------------------------------------- */ template <> inline void Material::unregisterInternal(InternalField & vect) { internal_vectors_real.erase(vect.getID()); } template <> inline void Material::unregisterInternal(InternalField & vect) { internal_vectors_uint.erase(vect.getID()); } template <> inline void Material::unregisterInternal(InternalField & vect) { internal_vectors_bool.erase(vect.getID()); } /* -------------------------------------------------------------------------- */ template inline bool Material::isInternal(__attribute__((unused)) const ID & id, __attribute__((unused)) const ElementKind & element_kind) const { AKANTU_TO_IMPLEMENT(); } template <> inline bool Material::isInternal(const ID & id, const ElementKind & element_kind) const { auto internal_array = internal_vectors_real.find(this->getID() + ":" + id); if (internal_array == internal_vectors_real.end() || internal_array->second->getElementKind() != element_kind) return false; return true; } /* -------------------------------------------------------------------------- */ template inline ElementTypeMap Material::getInternalDataPerElem(const ID & field_id, const ElementKind & element_kind) const { if (!this->template isInternal(field_id, element_kind)) AKANTU_EXCEPTION("Cannot find internal field " << id << " in material " << this->name); const InternalField & internal_field = this->template getInternal(field_id); const FEEngine & fe_engine = internal_field.getFEEngine(); UInt nb_data_per_quad = internal_field.getNbComponent(); ElementTypeMap res; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { using type_iterator = typename InternalField::type_iterator; type_iterator tit = internal_field.firstType(*gt); type_iterator tend = internal_field.lastType(*gt); for (; tit != tend; ++tit) { UInt nb_quadrature_points = fe_engine.getNbIntegrationPoints(*tit, *gt); res(*tit, *gt) = nb_data_per_quad * nb_quadrature_points; } } return res; } /* -------------------------------------------------------------------------- */ template void Material::flattenInternal(const std::string & field_id, ElementTypeMapArray & internal_flat, const GhostType ghost_type, ElementKind element_kind) const { if (!this->template isInternal(field_id, element_kind)) AKANTU_EXCEPTION("Cannot find internal field " << id << " in material " << this->name); const InternalField & internal_field = this->template getInternal(field_id); const FEEngine & fe_engine = internal_field.getFEEngine(); const Mesh & mesh = fe_engine.getMesh(); using type_iterator = typename InternalField::filter_type_iterator; type_iterator tit = internal_field.filterFirstType(ghost_type); type_iterator tend = internal_field.filterLastType(ghost_type); for (; tit != tend; ++tit) { ElementType type = *tit; const Array & src_vect = internal_field(type, ghost_type); const Array & filter = internal_field.getFilter(type, ghost_type); // total number of elements in the corresponding mesh UInt nb_element_dst = mesh.getNbElement(type, ghost_type); // number of element in the internal field UInt nb_element_src = filter.size(); // number of quadrature points per elem UInt nb_quad_per_elem = fe_engine.getNbIntegrationPoints(type); // number of data per quadrature point UInt nb_data_per_quad = internal_field.getNbComponent(); if (!internal_flat.exists(type, ghost_type)) { internal_flat.alloc(nb_element_dst * nb_quad_per_elem, nb_data_per_quad, type, ghost_type); } if (nb_element_src == 0) continue; // number of data per element UInt nb_data = nb_quad_per_elem * nb_data_per_quad; Array & dst_vect = internal_flat(type, ghost_type); dst_vect.resize(nb_element_dst * nb_quad_per_elem); Array::const_scalar_iterator it = filter.begin(); Array::const_scalar_iterator end = filter.end(); auto it_src = src_vect.begin_reinterpret(nb_data, nb_element_src); auto it_dst = dst_vect.begin_reinterpret(nb_data, nb_element_dst); for (; it != end; ++it, ++it_src) { it_dst[*it] = *it_src; } } } } // akantu #endif /* __AKANTU_MATERIAL_INLINE_IMPL_CC__ */ diff --git a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc index 15dcfb555..72e8c8fe5 100644 --- a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc +++ b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc @@ -1,266 +1,263 @@ /** * @file material_elastic_linear_anisotropic.cc * * @author Till Junge * @author Nicolas Richart * * @date creation: Wed Sep 25 2013 * @date last modification: Thu Oct 15 2015 * * @brief Anisotropic elastic material * * @section LICENSE * * Copyright (©) 2014, 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 . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "material_elastic_linear_anisotropic.hh" #include "solid_mechanics_model.hh" #include #include namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialElasticLinearAnisotropic::MaterialElasticLinearAnisotropic( SolidMechanicsModel & model, const ID & id, bool symmetric) : Material(model, id), rot_mat(dim, dim), Cprime(dim * dim, dim * dim), C(voigt_h::size, voigt_h::size), eigC(voigt_h::size), symmetric(symmetric), alpha(0), was_stiffness_assembled(false) { AKANTU_DEBUG_IN(); this->dir_vecs.push_back(std::make_unique>(dim)); (*this->dir_vecs.back())[0] = 1.; this->registerParam("n1", *(this->dir_vecs.back()), _pat_parsmod, "Direction of main material axis"); if (dim > 1) { this->dir_vecs.push_back(std::make_unique>(dim)); (*this->dir_vecs.back())[1] = 1.; this->registerParam("n2", *(this->dir_vecs.back()), _pat_parsmod, "Direction of secondary material axis"); } if (dim > 2) { this->dir_vecs.push_back(std::make_unique>(dim)); (*this->dir_vecs.back())[2] = 1.; this->registerParam("n3", *(this->dir_vecs.back()), _pat_parsmod, "Direction of tertiary material axis"); } for (UInt i = 0; i < voigt_h::size; ++i) { UInt start = 0; if (this->symmetric) { start = i; } for (UInt j = start; j < voigt_h::size; ++j) { std::stringstream param("C"); param << "C" << i + 1 << j + 1; this->registerParam(param.str(), this->Cprime(i, j), Real(0.), _pat_parsmod, "Coefficient " + param.str()); } } this->registerParam("alpha", this->alpha, _pat_parsmod, "Proportion of viscous stress"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); - - updateInternalParameters(); - AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::updateInternalParameters() { Material::updateInternalParameters(); if (this->symmetric) { for (UInt i = 0; i < voigt_h::size; ++i) { for (UInt j = i + 1; j < voigt_h::size; ++j) { this->Cprime(j, i) = this->Cprime(i, j); } } } this->rotateCprime(); this->C.eig(this->eigC); this->was_stiffness_assembled = false; } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::rotateCprime() { // start by filling the empty parts fo Cprime UInt diff = Dim * Dim - voigt_h::size; for (UInt i = voigt_h::size; i < Dim * Dim; ++i) { for (UInt j = 0; j < Dim * Dim; ++j) { this->Cprime(i, j) = this->Cprime(i - diff, j); } } for (UInt i = 0; i < Dim * Dim; ++i) { for (UInt j = voigt_h::size; j < Dim * Dim; ++j) { this->Cprime(i, j) = this->Cprime(i, j - diff); } } // construction of rotator tensor // normalise rotation matrix for (UInt j = 0; j < Dim; ++j) { Vector rot_vec = this->rot_mat(j); rot_vec = *this->dir_vecs[j]; rot_vec.normalize(); } // make sure the vectors form a right-handed base Vector test_axis(3); Vector v1(3), v2(3), v3(3); if (Dim == 2) { for (UInt i = 0; i < Dim; ++i) { v1[i] = this->rot_mat(0, i); v2[i] = this->rot_mat(1, i); v3[i] = 0.; } v3[2] = 1.; v1[2] = 0.; v2[2] = 0.; } else if (Dim == 3) { v1 = this->rot_mat(0); v2 = this->rot_mat(1); v3 = this->rot_mat(2); } test_axis.crossProduct(v1, v2); test_axis -= v3; if (test_axis.norm() > 8 * std::numeric_limits::epsilon()) { AKANTU_ERROR("The axis vectors do not form a right-handed coordinate " << "system. I. e., ||n1 x n2 - n3|| should be zero, but " << "it is " << test_axis.norm() << "."); } // create the rotator and the reverse rotator Matrix rotator(Dim * Dim, Dim * Dim); Matrix revrotor(Dim * Dim, Dim * Dim); for (UInt i = 0; i < Dim; ++i) { for (UInt j = 0; j < Dim; ++j) { for (UInt k = 0; k < Dim; ++k) { for (UInt l = 0; l < Dim; ++l) { UInt I = voigt_h::mat[i][j]; UInt J = voigt_h::mat[k][l]; rotator(I, J) = this->rot_mat(k, i) * this->rot_mat(l, j); revrotor(I, J) = this->rot_mat(i, k) * this->rot_mat(j, l); } } } } // create the full rotated matrix Matrix Cfull(Dim * Dim, Dim * Dim); Cfull = rotator * Cprime * revrotor; for (UInt i = 0; i < voigt_h::size; ++i) { for (UInt j = 0; j < voigt_h::size; ++j) { this->C(i, j) = Cfull(i, j); } } } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::computeStress(ElementType el_type, GhostType ghost_type) { // Wikipedia convention: // 2*eps_ij (i!=j) = voigt_eps_I // http://en.wikipedia.org/wiki/Voigt_notation AKANTU_DEBUG_IN(); Matrix strain(dim, dim); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); this->computeStressOnQuad(strain,sigma); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::computeTangentModuli( const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type) { AKANTU_DEBUG_IN(); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); this->computeTangentModuliOnQuad(tangent); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; this->was_stiffness_assembled = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::computePotentialEnergy( ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Material::computePotentialEnergy(el_type, ghost_type); AKANTU_DEBUG_ASSERT(!this->finite_deformation, "finite deformation not possible in material anisotropic " "(TO BE IMPLEMENTED)"); if (ghost_type != _not_ghost) return; Array::scalar_iterator epot = this->potential_energy(el_type, ghost_type).begin(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); computePotentialEnergyOnQuad(grad_u, sigma, *epot); ++epot; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Real MaterialElasticLinearAnisotropic::getCelerity( __attribute__((unused)) const Element & element) const { return std::sqrt(this->eigC(0) / rho); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(elastic_anisotropic, MaterialElasticLinearAnisotropic); } // akantu diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_damage_materials.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_damage_materials.cc index ed8a61267..cc5a0514b 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/test_damage_materials.cc +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_damage_materials.cc @@ -1,48 +1,48 @@ /* -------------------------------------------------------------------------- */ #include "material_marigo.hh" #include "material_mazars.hh" #include "solid_mechanics_model.hh" #include "test_material_fixtures.hh" #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; using types = ::testing::Types< Traits, Traits, Traits, Traits, Traits, Traits>; /*****************************************************************/ namespace { template class TestDamageMaterialFixture : public ::TestMaterialFixture {}; TYPED_TEST_CASE(TestDamageMaterialFixture, types); -TYPED_TEST(TestDamageMaterialFixture, DamageComputeStress) { +TYPED_TEST(TestDamageMaterialFixture, DISABLED_ComputeStress) { this->material->testComputeStress(); } -TYPED_TEST(TestDamageMaterialFixture, DamageEnergyDensity) { +TYPED_TEST(TestDamageMaterialFixture, DISABLED_EnergyDensity) { this->material->testEnergyDensity(); } -TYPED_TEST(TestDamageMaterialFixture, DamageComputeTangentModuli) { +TYPED_TEST(TestDamageMaterialFixture, DISABLED_ComputeTangentModuli) { this->material->testComputeTangentModuli(); } -TYPED_TEST(TestDamageMaterialFixture, DamageComputePushWaveSpeed) { +TYPED_TEST(TestDamageMaterialFixture, DISABLED_ComputePushWaveSpeed) { this->material->testPushWaveSpeed(); } -TYPED_TEST(TestDamageMaterialFixture, DamageComputeShearWaveSpeed) { +TYPED_TEST(TestDamageMaterialFixture, DISABLED_ComputeShearWaveSpeed) { this->material->testShearWaveSpeed(); } } /*****************************************************************/ diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_elastic_materials.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_elastic_materials.cc index ed3baeb2f..16e3330f6 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/test_elastic_materials.cc +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_elastic_materials.cc @@ -1,879 +1,807 @@ /* -------------------------------------------------------------------------- */ #include "material_elastic.hh" #include "material_elastic_orthotropic.hh" #include "solid_mechanics_model.hh" #include "test_material_fixtures.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; using types = ::testing::Types, Traits, Traits, Traits, Traits, Traits, Traits, Traits, Traits>; /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeStress() { Real E = 3.; setParam("E", E); Matrix eps = {{2}}; Matrix sigma(1, 1); Real sigma_th = 2; this->computeStressOnQuad(eps, sigma, sigma_th); auto solution = E * eps(0, 0) + sigma_th; ASSERT_NEAR(sigma(0, 0), solution, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testEnergyDensity() { Real eps = 2, sigma = 2; Real epot = 0; this->computePotentialEnergyOnQuad({{eps}}, {{sigma}}, epot); Real solution = 2; ASSERT_NEAR(epot, solution, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeTangentModuli() { Real E = 2; setParam("E", E); Matrix tangent(1, 1); this->computeTangentModuliOnQuad(tangent); ASSERT_NEAR(tangent(0, 0), E, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testPushWaveSpeed() { Real E = 3., rho = 2.; setParam("E", E); setParam("rho", rho); auto wave_speed = this->getPushWaveSpeed(Element()); auto solution = std::sqrt(E / rho); ASSERT_NEAR(wave_speed, solution, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testShearWaveSpeed() { ASSERT_THROW(this->getShearWaveSpeed(Element()), debug::Exception); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeStress() { Real E = 1.; Real nu = .3; Real sigma_th = 0.3; // thermal stress setParam("E", E); setParam("nu", nu); Matrix rotation_matrix = getRandomRotation2d(); auto grad_u = this->getDeviatoricStrain(1.).block(0, 0, 2, 2); auto grad_u_rot = this->applyRotation(grad_u, rotation_matrix); Matrix sigma_rot(2, 2); this->computeStressOnQuad(grad_u_rot, sigma_rot, sigma_th); auto sigma = this->reverseRotation(sigma_rot, rotation_matrix); Matrix identity(2, 2); identity.eye(); Matrix sigma_expected = - 0.5 * E / (1 + nu) * (grad_u + grad_u.transpose()) + sigma_th * identity; + 0.5 * E / (1 + nu) * (grad_u + grad_u.transpose()) + sigma_th * identity; auto diff = sigma - sigma_expected; Real stress_error = diff.norm(); ASSERT_NEAR(stress_error, 0., 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testEnergyDensity() { Matrix sigma = {{1, 2}, {2, 4}}; Matrix eps = {{1, 0}, {0, 1}}; Real epot = 0; Real solution = 2.5; this->computePotentialEnergyOnQuad(eps, sigma, epot); ASSERT_NEAR(epot, solution, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeTangentModuli() { Real E = 1.; Real nu = .3; setParam("E", E); setParam("nu", nu); Matrix tangent(3, 3); /* Plane Strain */ // clang-format off Matrix solution = { {1 - nu, nu, 0}, {nu, 1 - nu, 0}, {0, 0, (1 - 2 * nu) / 2}, }; // clang-format on solution *= E / ((1 + nu) * (1 - 2 * nu)); this->computeTangentModuliOnQuad(tangent); Real tangent_error = (tangent - solution).norm(); ASSERT_NEAR(tangent_error, 0, 1e-14); /* Plane Stress */ this->plane_stress = true; this->updateInternalParameters(); // clang-format off solution = { {1, nu, 0}, {nu, 1, 0}, {0, 0, (1 - nu) / 2}, }; // clang-format on solution *= E / (1 - nu * nu); this->computeTangentModuliOnQuad(tangent); tangent_error = (tangent - solution).norm(); ASSERT_NEAR(tangent_error, 0, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testPushWaveSpeed() { Real E = 1.; Real nu = .3; Real rho = 2; setParam("E", E); setParam("nu", nu); setParam("rho", rho); auto wave_speed = this->getPushWaveSpeed(Element()); Real K = E / (3 * (1 - 2 * nu)); Real mu = E / (2 * (1 + nu)); Real sol = std::sqrt((K + 4. / 3 * mu) / rho); ASSERT_NEAR(wave_speed, sol, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testShearWaveSpeed() { Real E = 1.; Real nu = .3; Real rho = 2; setParam("E", E); setParam("nu", nu); setParam("rho", rho); auto wave_speed = this->getShearWaveSpeed(Element()); Real mu = E / (2 * (1 + nu)); Real sol = std::sqrt(mu / rho); ASSERT_NEAR(wave_speed, sol, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeStress() { Real E = 1.; Real nu = .3; Real sigma_th = 0.3; // thermal stress setParam("E", E); setParam("nu", nu); Matrix rotation_matrix = getRandomRotation3d(); auto grad_u = this->getDeviatoricStrain(1.); auto grad_u_rot = this->applyRotation(grad_u, rotation_matrix); Matrix sigma_rot(3, 3); this->computeStressOnQuad(grad_u_rot, sigma_rot, sigma_th); auto sigma = this->reverseRotation(sigma_rot, rotation_matrix); Matrix identity(3, 3); identity.eye(); Matrix sigma_expected = 0.5 * E / (1 + nu) * (grad_u + grad_u.transpose()) + sigma_th * identity; auto diff = sigma - sigma_expected; Real stress_error = diff.norm(); ASSERT_NEAR(stress_error, 0., 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testEnergyDensity() { Matrix sigma = {{1, 2, 3}, {2, 4, 5}, {3, 5, 6}}; Matrix eps = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}; Real epot = 0; Real solution = 5.5; this->computePotentialEnergyOnQuad(eps, sigma, epot); ASSERT_NEAR(epot, solution, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeTangentModuli() { Real E = 1.; Real nu = .3; setParam("E", E); setParam("nu", nu); Matrix tangent(6, 6); // clang-format off Matrix solution = { {1 - nu, nu, nu, 0, 0, 0}, {nu, 1 - nu, nu, 0, 0, 0}, {nu, nu, 1 - nu, 0, 0, 0}, {0, 0, 0, (1 - 2 * nu) / 2, 0, 0}, {0, 0, 0, 0, (1 - 2 * nu) / 2, 0}, {0, 0, 0, 0, 0, (1 - 2 * nu) / 2}, }; // clang-format on solution *= E / ((1 + nu) * (1 - 2 * nu)); this->computeTangentModuliOnQuad(tangent); Real tangent_error = (tangent - solution).norm(); ASSERT_NEAR(tangent_error, 0, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testPushWaveSpeed() { Real E = 1.; Real nu = .3; Real rho = 2; setParam("E", E); setParam("nu", nu); setParam("rho", rho); auto wave_speed = this->getPushWaveSpeed(Element()); Real K = E / (3 * (1 - 2 * nu)); Real mu = E / (2 * (1 + nu)); Real sol = std::sqrt((K + 4. / 3 * mu) / rho); ASSERT_NEAR(wave_speed, sol, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testShearWaveSpeed() { Real E = 1.; Real nu = .3; Real rho = 2; setParam("E", E); setParam("nu", nu); setParam("rho", rho); auto wave_speed = this->getShearWaveSpeed(Element()); Real mu = E / (2 * (1 + nu)); Real sol = std::sqrt(mu / rho); ASSERT_NEAR(wave_speed, sol, 1e-14); } /* -------------------------------------------------------------------------- */ -template <> void FriendMaterial>::testComputeStress() { +template <> +void FriendMaterial>::testComputeStress() { UInt Dim = 2; - Real E1 = 1.; - Real E2 = 2.; - Real nu12 = 0.1; - Real G12 = 2.; - setParamNoUpdate("E1", E1); - setParamNoUpdate("E2", E2); - setParamNoUpdate("nu12", nu12); - setParamNoUpdate("G12", G12); + this->E1 = 1.; + this->E2 = 2.; + this->nu12 = 0.1; + this->G12 = 2.; // material frame of reference is rotate by rotation_matrix starting from // canonical basis Matrix rotation_matrix = getRandomRotation2d(); // canonical basis as expressed in the material frame of reference, as // required by MaterialElasticOrthotropic class (it is simply given by the // columns of the rotation_matrix; the lines give the material basis expressed // in the canonical frame of reference) - Vector v1(Dim); - Vector v2(Dim); - for (UInt i = 0; i < Dim; ++i) { - v1[i] = rotation_matrix(i, 0); - v2[i] = rotation_matrix(i, 1); - } - - setParamNoUpdate("n1", v1.normalize()); - setParamNoUpdate("n2", v2.normalize()); + *this->dir_vecs[0] = rotation_matrix(0); + *this->dir_vecs[1] = rotation_matrix(1); // set internal Cijkl matrix expressed in the canonical frame of reference this->updateInternalParameters(); // gradient in material frame of reference - auto grad_u = ( this->getDeviatoricStrain(2.) + this->getHydrostaticStrain(1.) ).block(0, 0, 2, 2); + auto grad_u = (this->getDeviatoricStrain(2.) + this->getHydrostaticStrain(1.)) + .block(0, 0, 2, 2); -// gradient in canonical basis (we need to rotate *back* to the canonical -// basis) + // gradient in canonical basis (we need to rotate *back* to the canonical + // basis) auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix); -// stress in the canonical basis + // stress in the canonical basis Matrix sigma_rot(2, 2); this->computeStressOnQuad(grad_u_rot, sigma_rot); -// stress in the material reference (we need to apply the rotation) + // stress in the material reference (we need to apply the rotation) auto sigma = this->applyRotation(sigma_rot, rotation_matrix); -// construction of Cijkl engineering tensor in the *material* frame of reference -// ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13 - Real nu21 = nu12*E2/E1; - Real gamma = 1 / ( 1 - nu12*nu21 ); + // construction of Cijkl engineering tensor in the *material* frame of + // reference + // ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13 + Real nu21 = nu12 * E2 / E1; + Real gamma = 1 / (1 - nu12 * nu21); - Matrix C_expected(2*Dim,2*Dim,0); - C_expected(0,0) = gamma * E1 ; - C_expected(1,1) = gamma * E2 ; - C_expected(2,2) = G12; + Matrix C_expected(2 * Dim, 2 * Dim, 0); + C_expected(0, 0) = gamma * E1; + C_expected(1, 1) = gamma * E2; + C_expected(2, 2) = G12; - C_expected(1,0) = C_expected(0,1) = gamma * E1 * nu21; + C_expected(1, 0) = C_expected(0, 1) = gamma * E1 * nu21; -// epsilon is computed directly in the *material* frame of reference - Matrix epsilon = 0.5 * ( grad_u + grad_u.transpose() ); + // epsilon is computed directly in the *material* frame of reference + Matrix epsilon = 0.5 * (grad_u + grad_u.transpose()); -// sigma_expected is computed directly in the *material* frame of reference + // sigma_expected is computed directly in the *material* frame of reference Matrix sigma_expected(Dim, Dim); - for ( UInt i = 0; i < Dim; ++i ) { - for ( UInt j = 0; j < Dim; ++j ) { - sigma_expected(i,i) += C_expected(i,j)*epsilon(j,j); + for (UInt i = 0; i < Dim; ++i) { + for (UInt j = 0; j < Dim; ++j) { + sigma_expected(i, i) += C_expected(i, j) * epsilon(j, j); } } - sigma_expected(0,1) = sigma_expected(1,0) = C_expected(2,2) * 2 * epsilon(0,1); + sigma_expected(0, 1) = sigma_expected(1, 0) = + C_expected(2, 2) * 2 * epsilon(0, 1); // sigmas are checked in the *material* frame of reference auto diff = sigma - sigma_expected; Real stress_error = diff.norm(); ASSERT_NEAR(stress_error, 0., 1e-13); } /* -------------------------------------------------------------------------- */ -template <> void FriendMaterial>::testEnergyDensity() { +template <> +void FriendMaterial>::testEnergyDensity() { Matrix sigma = {{1, 2}, {2, 4}}; Matrix eps = {{1, 0}, {0, 1}}; Real epot = 0; Real solution = 2.5; this->computePotentialEnergyOnQuad(eps, sigma, epot); ASSERT_NEAR(epot, solution, 1e-14); } /* -------------------------------------------------------------------------- */ -template <> void FriendMaterial>::testComputeTangentModuli() { +template <> +void FriendMaterial>::testComputeTangentModuli() { // Note: for this test material and canonical basis coincide Vector n1 = {1, 0}; Vector n2 = {0, 1}; Real E1 = 1.; Real E2 = 2.; Real nu12 = 0.1; Real G12 = 2.; - setParamNoUpdate("n1", n1); - setParamNoUpdate("n2", n2); - setParamNoUpdate("E1", E1); - setParamNoUpdate("E2", E2); - setParamNoUpdate("nu12", nu12); - setParamNoUpdate("G12", G12); + *this->dir_vecs[0] = n1; + *this->dir_vecs[1] = n2; + + this->E1 = E1; + this->E2 = E2; + this->nu12 = nu12; + this->G12 = G12; // set internal Cijkl matrix expressed in the canonical frame of reference this->updateInternalParameters(); -// construction of Cijkl engineering tensor in the *material* frame of reference -// ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13 - Real nu21 = nu12*E2/E1; - Real gamma = 1 / ( 1 - nu12*nu21 ); + // construction of Cijkl engineering tensor in the *material* frame of + // reference + // ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13 + Real nu21 = nu12 * E2 / E1; + Real gamma = 1 / (1 - nu12 * nu21); - Matrix C_expected(3,3); - C_expected(0,0) = gamma * E1 ; - C_expected(1,1) = gamma * E2 ; - C_expected(2,2) = G12; + Matrix C_expected(3, 3); + C_expected(0, 0) = gamma * E1; + C_expected(1, 1) = gamma * E2; + C_expected(2, 2) = G12; - C_expected(1,0) = C_expected(0,1) = gamma * E1 * nu21; + C_expected(1, 0) = C_expected(0, 1) = gamma * E1 * nu21; - Matrix tangent(3,3); + Matrix tangent(3, 3); this->computeTangentModuliOnQuad(tangent); Real tangent_error = (tangent - C_expected).norm(); ASSERT_NEAR(tangent_error, 0, 1e-14); } /* -------------------------------------------------------------------------- */ -template <> void FriendMaterial>::testComputeStress() { +template <> +void FriendMaterial>::testComputeStress() { UInt Dim = 3; Real E1 = 1.; Real E2 = 2.; Real E3 = 3.; Real nu12 = 0.1; Real nu13 = 0.2; Real nu23 = 0.3; Real G12 = 2.; Real G13 = 3.; Real G23 = 1.; - setParamNoUpdate("E1", E1); - setParamNoUpdate("E2", E2); - setParamNoUpdate("E3", E3); - setParamNoUpdate("nu12", nu12); - setParamNoUpdate("nu13", nu13); - setParamNoUpdate("nu23", nu23); - setParamNoUpdate("G12", G12); - setParamNoUpdate("G13", G13); - setParamNoUpdate("G23", G23); + this->E1 = E1; + this->E2 = E2; + this->E3 = E3; + this->nu12 = nu12; + this->nu13 = nu13; + this->nu23 = nu23; + this->G12 = G12; + this->G13 = G13; + this->G23 = G23; // material frame of reference is rotate by rotation_matrix starting from // canonical basis Matrix rotation_matrix = getRandomRotation3d(); // canonical basis as expressed in the material frame of reference, as // required by MaterialElasticOrthotropic class (it is simply given by the // columns of the rotation_matrix; the lines give the material basis expressed // in the canonical frame of reference) - Vector v1(Dim); - Vector v2(Dim); - Vector v3(Dim); - for (UInt i = 0; i < Dim; ++i) { - v1[i] = rotation_matrix(i, 0); - v2[i] = rotation_matrix(i, 1); - v3[i] = rotation_matrix(i, 2); - } - - setParamNoUpdate("n1", v1.normalize()); - setParamNoUpdate("n2", v2.normalize()); - setParamNoUpdate("n3", v3.normalize()); + *this->dir_vecs[0] = rotation_matrix(0); + *this->dir_vecs[1] = rotation_matrix(1); + *this->dir_vecs[2] = rotation_matrix(2); // set internal Cijkl matrix expressed in the canonical frame of reference this->updateInternalParameters(); // gradient in material frame of reference auto grad_u = this->getDeviatoricStrain(2.) + this->getHydrostaticStrain(1.); // gradient in canonical basis (we need to rotate *back* to the canonical // basis) auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix); // stress in the canonical basis Matrix sigma_rot(3, 3); this->computeStressOnQuad(grad_u_rot, sigma_rot); // stress in the material reference (we need to apply the rotation) auto sigma = this->applyRotation(sigma_rot, rotation_matrix); -// construction of Cijkl engineering tensor in the *material* frame of reference -// ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13 - Real nu21 = nu12*E2/E1; - Real nu31 = nu13*E3/E1; - Real nu32 = nu23*E3/E2; - Real gamma = 1 / ( 1 - nu12*nu21 - nu23*nu32 - nu31*nu13 - 2*nu21*nu32*nu13 ); + // construction of Cijkl engineering tensor in the *material* frame of + // reference + // ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13 + Real nu21 = nu12 * E2 / E1; + Real nu31 = nu13 * E3 / E1; + Real nu32 = nu23 * E3 / E2; + Real gamma = 1 / (1 - nu12 * nu21 - nu23 * nu32 - nu31 * nu13 - + 2 * nu21 * nu32 * nu13); - Matrix C_expected(6,6); - C_expected(0,0) = gamma * E1 * ( 1-nu23*nu32 ); - C_expected(1,1) = gamma * E2 * ( 1-nu13*nu31 ); - C_expected(2,2) = gamma * E3 * ( 1-nu12*nu21 ); + Matrix C_expected(6, 6); + C_expected(0, 0) = gamma * E1 * (1 - nu23 * nu32); + C_expected(1, 1) = gamma * E2 * (1 - nu13 * nu31); + C_expected(2, 2) = gamma * E3 * (1 - nu12 * nu21); - C_expected(1,0) = C_expected(0,1) = gamma * E1 * ( nu21 + nu31*nu23 ); - C_expected(2,0) = C_expected(0,2) = gamma * E1 * ( nu31 + nu21*nu32 ); - C_expected(2,1) = C_expected(1,2) = gamma * E2 * ( nu32 + nu12*nu31 ); + C_expected(1, 0) = C_expected(0, 1) = gamma * E1 * (nu21 + nu31 * nu23); + C_expected(2, 0) = C_expected(0, 2) = gamma * E1 * (nu31 + nu21 * nu32); + C_expected(2, 1) = C_expected(1, 2) = gamma * E2 * (nu32 + nu12 * nu31); - C_expected(3,3) = G23; - C_expected(4,4) = G13; - C_expected(5,5) = G12; + C_expected(3, 3) = G23; + C_expected(4, 4) = G13; + C_expected(5, 5) = G12; // epsilon is computed directly in the *material* frame of reference - Matrix epsilon = 0.5 * ( grad_u + grad_u.transpose() ); + Matrix epsilon = 0.5 * (grad_u + grad_u.transpose()); // sigma_expected is computed directly in the *material* frame of reference Matrix sigma_expected(Dim, Dim); - for ( UInt i = 0; i < Dim; ++i ) { - for ( UInt j = 0; j < Dim; ++j ) { - sigma_expected(i,i) += C_expected(i,j)*epsilon(j,j); + for (UInt i = 0; i < Dim; ++i) { + for (UInt j = 0; j < Dim; ++j) { + sigma_expected(i, i) += C_expected(i, j) * epsilon(j, j); } } - sigma_expected(0,1) = C_expected(5,5) * 2 * epsilon(0,1); - sigma_expected(0,2) = C_expected(4,4) * 2 * epsilon(0,2); - sigma_expected(1,2) = C_expected(3,3) * 2 * epsilon(1,2); - sigma_expected(1,0) = sigma_expected(0,1); - sigma_expected(2,0) = sigma_expected(0,2); - sigma_expected(2,1) = sigma_expected(1,2); + sigma_expected(0, 1) = C_expected(5, 5) * 2 * epsilon(0, 1); + sigma_expected(0, 2) = C_expected(4, 4) * 2 * epsilon(0, 2); + sigma_expected(1, 2) = C_expected(3, 3) * 2 * epsilon(1, 2); + sigma_expected(1, 0) = sigma_expected(0, 1); + sigma_expected(2, 0) = sigma_expected(0, 2); + sigma_expected(2, 1) = sigma_expected(1, 2); // sigmas are checked in the *material* frame of reference auto diff = sigma - sigma_expected; Real stress_error = diff.norm(); ASSERT_NEAR(stress_error, 0., 1e-13); } /* -------------------------------------------------------------------------- */ -template <> void FriendMaterial>::testEnergyDensity() { +template <> +void FriendMaterial>::testEnergyDensity() { Matrix sigma = {{1, 2, 3}, {2, 4, 5}, {3, 5, 6}}; Matrix eps = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}; Real epot = 0; Real solution = 5.5; this->computePotentialEnergyOnQuad(eps, sigma, epot); ASSERT_NEAR(epot, solution, 1e-14); } /* -------------------------------------------------------------------------- */ -template <> void FriendMaterial>::testComputeTangentModuli() { +template <> +void FriendMaterial>::testComputeTangentModuli() { // Note: for this test material and canonical basis coincide UInt Dim = 3; Vector n1 = {1, 0, 0}; Vector n2 = {0, 1, 0}; Vector n3 = {0, 0, 1}; Real E1 = 1.; Real E2 = 2.; Real E3 = 3.; Real nu12 = 0.1; Real nu13 = 0.2; Real nu23 = 0.3; Real G12 = 2.; Real G13 = 3.; Real G23 = 1.; - setParamNoUpdate("n1", n1); - setParamNoUpdate("n2", n2); - setParamNoUpdate("n3", n3); - setParamNoUpdate("E1", E1); - setParamNoUpdate("E2", E2); - setParamNoUpdate("E3", E3); - setParamNoUpdate("nu12", nu12); - setParamNoUpdate("nu13", nu13); - setParamNoUpdate("nu23", nu23); - setParamNoUpdate("G12", G12); - setParamNoUpdate("G13", G13); - setParamNoUpdate("G23", G23); + *this->dir_vecs[0] = n1; + *this->dir_vecs[1] = n2; + *this->dir_vecs[2] = n3; + this->E1 = E1; + this->E2 = E2; + this->E3 = E3; + this->nu12 = nu12; + this->nu13 = nu13; + this->nu23 = nu23; + this->G12 = G12; + this->G13 = G13; + this->G23 = G23; // set internal Cijkl matrix expressed in the canonical frame of reference this->updateInternalParameters(); -// construction of Cijkl engineering tensor in the *material* frame of reference -// ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13 - Real nu21 = nu12*E2/E1; - Real nu31 = nu13*E3/E1; - Real nu32 = nu23*E3/E2; - Real gamma = 1 / ( 1 - nu12*nu21 - nu23*nu32 - nu31*nu13 - 2*nu21*nu32*nu13 ); + // construction of Cijkl engineering tensor in the *material* frame of + // reference + // ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13 + Real nu21 = nu12 * E2 / E1; + Real nu31 = nu13 * E3 / E1; + Real nu32 = nu23 * E3 / E2; + Real gamma = 1 / (1 - nu12 * nu21 - nu23 * nu32 - nu31 * nu13 - + 2 * nu21 * nu32 * nu13); - Matrix C_expected(2*Dim,2*Dim,0); - C_expected(0,0) = gamma * E1 * ( 1-nu23*nu32 ); - C_expected(1,1) = gamma * E2 * ( 1-nu13*nu31 ); - C_expected(2,2) = gamma * E3 * ( 1-nu12*nu21 ); + Matrix C_expected(2 * Dim, 2 * Dim, 0); + C_expected(0, 0) = gamma * E1 * (1 - nu23 * nu32); + C_expected(1, 1) = gamma * E2 * (1 - nu13 * nu31); + C_expected(2, 2) = gamma * E3 * (1 - nu12 * nu21); - C_expected(1,0) = C_expected(0,1) = gamma * E1 * ( nu21 + nu31*nu23 ); - C_expected(2,0) = C_expected(0,2) = gamma * E1 * ( nu31 + nu21*nu32 ); - C_expected(2,1) = C_expected(1,2) = gamma * E2 * ( nu32 + nu12*nu31 ); + C_expected(1, 0) = C_expected(0, 1) = gamma * E1 * (nu21 + nu31 * nu23); + C_expected(2, 0) = C_expected(0, 2) = gamma * E1 * (nu31 + nu21 * nu32); + C_expected(2, 1) = C_expected(1, 2) = gamma * E2 * (nu32 + nu12 * nu31); - C_expected(3,3) = G23; - C_expected(4,4) = G13; - C_expected(5,5) = G12; + C_expected(3, 3) = G23; + C_expected(4, 4) = G13; + C_expected(5, 5) = G12; - Matrix tangent(6,6); + Matrix tangent(6, 6); this->computeTangentModuliOnQuad(tangent); Real tangent_error = (tangent - C_expected).norm(); ASSERT_NEAR(tangent_error, 0, 1e-14); } /* -------------------------------------------------------------------------- */ -template <> void FriendMaterial>::testComputeStress() { +template <> +void FriendMaterial>::testComputeStress() { UInt Dim = 2; Matrix C = { - {1.0, 0.3, 0.4}, - {0.3, 2.0, 0.1}, - {0.4, 0.1, 1.5}, + {1.0, 0.3, 0.4}, {0.3, 2.0, 0.1}, {0.4, 0.1, 1.5}, }; - setParamNoUpdate("C11", C(0,0)); - setParamNoUpdate("C12", C(0,1)); - setParamNoUpdate("C13", C(0,2)); - setParamNoUpdate("C22", C(1,1)); - setParamNoUpdate("C23", C(1,2)); - setParamNoUpdate("C33", C(2,2)); + this->Cprime = C; // material frame of reference is rotate by rotation_matrix starting from // canonical basis Matrix rotation_matrix = getRandomRotation2d(); // canonical basis as expressed in the material frame of reference, as // required by MaterialElasticLinearAnisotropic class (it is simply given by // the columns of the rotation_matrix; the lines give the material basis // expressed in the canonical frame of reference) - Vector v1(Dim); - Vector v2(Dim); - for (UInt i = 0; i < Dim; ++i) { - v1[i] = rotation_matrix(i, 0); - v2[i] = rotation_matrix(i, 1); - } - - setParamNoUpdate("n1", v1.normalize()); - setParamNoUpdate("n2", v2.normalize()); + *this->dir_vecs[0] = rotation_matrix(0); + *this->dir_vecs[1] = rotation_matrix(1); // set internal Cijkl matrix expressed in the canonical frame of reference this->updateInternalParameters(); // gradient in material frame of reference - auto grad_u = (this->getDeviatoricStrain(2.) + this->getHydrostaticStrain(1.)).block(0, 0, 2, 2); + auto grad_u = (this->getDeviatoricStrain(2.) + this->getHydrostaticStrain(1.)) + .block(0, 0, 2, 2); // gradient in canonical basis (we need to rotate *back* to the canonical // basis) auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix); // stress in the canonical basis Matrix sigma_rot(2, 2); this->computeStressOnQuad(grad_u_rot, sigma_rot); // stress in the material reference (we need to apply the rotation) auto sigma = this->applyRotation(sigma_rot, rotation_matrix); // epsilon is computed directly in the *material* frame of reference - Matrix epsilon = 0.5 * ( grad_u + grad_u.transpose() ); + Matrix epsilon = 0.5 * (grad_u + grad_u.transpose()); Vector epsilon_voigt(3); - epsilon_voigt(0) = epsilon(0,0); - epsilon_voigt(1) = epsilon(1,1); - epsilon_voigt(2) = 2 * epsilon(0,1); + epsilon_voigt(0) = epsilon(0, 0); + epsilon_voigt(1) = epsilon(1, 1); + epsilon_voigt(2) = 2 * epsilon(0, 1); // sigma_expected is computed directly in the *material* frame of reference Vector sigma_voigt = C * epsilon_voigt; Matrix sigma_expected(Dim, Dim); - sigma_expected(0,0) = sigma_voigt(0); - sigma_expected(1,1) = sigma_voigt(1); - sigma_expected(0,1) = sigma_expected(1,0) = sigma_voigt(2); + sigma_expected(0, 0) = sigma_voigt(0); + sigma_expected(1, 1) = sigma_voigt(1); + sigma_expected(0, 1) = sigma_expected(1, 0) = sigma_voigt(2); // sigmas are checked in the *material* frame of reference auto diff = sigma - sigma_expected; Real stress_error = diff.norm(); ASSERT_NEAR(stress_error, 0., 1e-13); } /* -------------------------------------------------------------------------- */ -template <> void FriendMaterial>::testEnergyDensity() { +template <> +void FriendMaterial>::testEnergyDensity() { Matrix sigma = {{1, 2}, {2, 4}}; Matrix eps = {{1, 0}, {0, 1}}; Real epot = 0; Real solution = 2.5; this->computePotentialEnergyOnQuad(eps, sigma, epot); ASSERT_NEAR(epot, solution, 1e-14); } /* -------------------------------------------------------------------------- */ -template <> void FriendMaterial>::testComputeTangentModuli() { +template <> +void FriendMaterial< + MaterialElasticLinearAnisotropic<2>>::testComputeTangentModuli() { // Note: for this test material and canonical basis coincide Matrix C = { - {1.0, 0.3, 0.4}, - {0.3, 2.0, 0.1}, - {0.4, 0.1, 1.5}, + {1.0, 0.3, 0.4}, {0.3, 2.0, 0.1}, {0.4, 0.1, 1.5}, }; - setParamNoUpdate("C11", C(0,0)); - setParamNoUpdate("C12", C(0,1)); - setParamNoUpdate("C13", C(0,2)); - setParamNoUpdate("C22", C(1,1)); - setParamNoUpdate("C23", C(1,2)); - setParamNoUpdate("C33", C(2,2)); + this->Cprime = C; // set internal Cijkl matrix expressed in the canonical frame of reference this->updateInternalParameters(); - Matrix tangent(3,3); + Matrix tangent(3, 3); this->computeTangentModuliOnQuad(tangent); Real tangent_error = (tangent - C).norm(); ASSERT_NEAR(tangent_error, 0, 1e-14); } /* -------------------------------------------------------------------------- */ -template <> void FriendMaterial>::testComputeStress() { +template <> +void FriendMaterial>::testComputeStress() { UInt Dim = 3; Matrix C = { - {1.0, 0.3, 0.4, 0.3, 0.2, 0.1}, - {0.3, 2.0, 0.1, 0.2, 0.3, 0.2}, - {0.4, 0.1, 1.5, 0.1, 0.4, 0.3}, - {0.3, 0.2, 0.1, 2.4, 0.1, 0.4}, - {0.2, 0.3, 0.4, 0.1, 0.9, 0.1}, - {0.1, 0.2, 0.3, 0.4, 0.1, 1.2}, + {1.0, 0.3, 0.4, 0.3, 0.2, 0.1}, {0.3, 2.0, 0.1, 0.2, 0.3, 0.2}, + {0.4, 0.1, 1.5, 0.1, 0.4, 0.3}, {0.3, 0.2, 0.1, 2.4, 0.1, 0.4}, + {0.2, 0.3, 0.4, 0.1, 0.9, 0.1}, {0.1, 0.2, 0.3, 0.4, 0.1, 1.2}, }; - setParamNoUpdate("C11", C(0,0)); - setParamNoUpdate("C12", C(0,1)); - setParamNoUpdate("C13", C(0,2)); - setParamNoUpdate("C14", C(0,3)); - setParamNoUpdate("C15", C(0,4)); - setParamNoUpdate("C16", C(0,5)); - setParamNoUpdate("C22", C(1,1)); - setParamNoUpdate("C23", C(1,2)); - setParamNoUpdate("C24", C(1,3)); - setParamNoUpdate("C25", C(1,4)); - setParamNoUpdate("C26", C(1,5)); - setParamNoUpdate("C33", C(2,2)); - setParamNoUpdate("C34", C(2,3)); - setParamNoUpdate("C35", C(2,4)); - setParamNoUpdate("C36", C(2,5)); - setParamNoUpdate("C44", C(3,3)); - setParamNoUpdate("C45", C(3,4)); - setParamNoUpdate("C46", C(3,5)); - setParamNoUpdate("C55", C(4,4)); - setParamNoUpdate("C56", C(4,5)); - setParamNoUpdate("C66", C(5,5)); + this->Cprime = C; // material frame of reference is rotate by rotation_matrix starting from // canonical basis Matrix rotation_matrix = getRandomRotation3d(); // canonical basis as expressed in the material frame of reference, as // required by MaterialElasticLinearAnisotropic class (it is simply given by // the columns of the rotation_matrix; the lines give the material basis // expressed in the canonical frame of reference) - Vector v1(Dim); - Vector v2(Dim); - Vector v3(Dim); - for (UInt i = 0; i < Dim; ++i) { - v1[i] = rotation_matrix(i, 0); - v2[i] = rotation_matrix(i, 1); - v3[i] = rotation_matrix(i, 2); - } - - setParamNoUpdate("n1", v1.normalize()); - setParamNoUpdate("n2", v2.normalize()); - setParamNoUpdate("n3", v3.normalize()); + *this->dir_vecs[0] = rotation_matrix(0); + *this->dir_vecs[1] = rotation_matrix(1); + *this->dir_vecs[2] = rotation_matrix(2); // set internal Cijkl matrix expressed in the canonical frame of reference this->updateInternalParameters(); // gradient in material frame of reference auto grad_u = this->getDeviatoricStrain(2.) + this->getHydrostaticStrain(1.); // gradient in canonical basis (we need to rotate *back* to the canonical // basis) auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix); // stress in the canonical basis Matrix sigma_rot(3, 3); this->computeStressOnQuad(grad_u_rot, sigma_rot); // stress in the material reference (we need to apply the rotation) auto sigma = this->applyRotation(sigma_rot, rotation_matrix); // epsilon is computed directly in the *material* frame of reference - Matrix epsilon = 0.5 * ( grad_u + grad_u.transpose() ); + Matrix epsilon = 0.5 * (grad_u + grad_u.transpose()); Vector epsilon_voigt(6); - epsilon_voigt(0) = epsilon(0,0); - epsilon_voigt(1) = epsilon(1,1); - epsilon_voigt(2) = epsilon(2,2); - epsilon_voigt(3) = 2 * epsilon(1,2); - epsilon_voigt(4) = 2 * epsilon(0,2); - epsilon_voigt(5) = 2 * epsilon(0,1); + epsilon_voigt(0) = epsilon(0, 0); + epsilon_voigt(1) = epsilon(1, 1); + epsilon_voigt(2) = epsilon(2, 2); + epsilon_voigt(3) = 2 * epsilon(1, 2); + epsilon_voigt(4) = 2 * epsilon(0, 2); + epsilon_voigt(5) = 2 * epsilon(0, 1); // sigma_expected is computed directly in the *material* frame of reference Vector sigma_voigt = C * epsilon_voigt; Matrix sigma_expected(Dim, Dim); - sigma_expected(0,0) = sigma_voigt(0); - sigma_expected(1,1) = sigma_voigt(1); - sigma_expected(2,2) = sigma_voigt(2); - sigma_expected(1,2) = sigma_expected(2,1) = sigma_voigt(3); - sigma_expected(0,2) = sigma_expected(2,0) = sigma_voigt(4); - sigma_expected(0,1) = sigma_expected(1,0) = sigma_voigt(5); + sigma_expected(0, 0) = sigma_voigt(0); + sigma_expected(1, 1) = sigma_voigt(1); + sigma_expected(2, 2) = sigma_voigt(2); + sigma_expected(1, 2) = sigma_expected(2, 1) = sigma_voigt(3); + sigma_expected(0, 2) = sigma_expected(2, 0) = sigma_voigt(4); + sigma_expected(0, 1) = sigma_expected(1, 0) = sigma_voigt(5); // sigmas are checked in the *material* frame of reference auto diff = sigma - sigma_expected; Real stress_error = diff.norm(); ASSERT_NEAR(stress_error, 0., 1e-13); } /* -------------------------------------------------------------------------- */ -template <> void FriendMaterial>::testEnergyDensity() { +template <> +void FriendMaterial>::testEnergyDensity() { Matrix sigma = {{1, 2, 3}, {2, 4, 5}, {3, 5, 6}}; Matrix eps = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}; Real epot = 0; Real solution = 5.5; this->computePotentialEnergyOnQuad(eps, sigma, epot); ASSERT_NEAR(epot, solution, 1e-14); } /* -------------------------------------------------------------------------- */ -template <> void FriendMaterial>::testComputeTangentModuli() { +template <> +void FriendMaterial< + MaterialElasticLinearAnisotropic<3>>::testComputeTangentModuli() { // Note: for this test material and canonical basis coincide Matrix C = { - {1.0, 0.3, 0.4, 0.3, 0.2, 0.1}, - {0.3, 2.0, 0.1, 0.2, 0.3, 0.2}, - {0.4, 0.1, 1.5, 0.1, 0.4, 0.3}, - {0.3, 0.2, 0.1, 2.4, 0.1, 0.4}, - {0.2, 0.3, 0.4, 0.1, 0.9, 0.1}, - {0.1, 0.2, 0.3, 0.4, 0.1, 1.2}, + {1.0, 0.3, 0.4, 0.3, 0.2, 0.1}, {0.3, 2.0, 0.1, 0.2, 0.3, 0.2}, + {0.4, 0.1, 1.5, 0.1, 0.4, 0.3}, {0.3, 0.2, 0.1, 2.4, 0.1, 0.4}, + {0.2, 0.3, 0.4, 0.1, 0.9, 0.1}, {0.1, 0.2, 0.3, 0.4, 0.1, 1.2}, }; - setParamNoUpdate("C11", C(0,0)); - setParamNoUpdate("C12", C(0,1)); - setParamNoUpdate("C13", C(0,2)); - setParamNoUpdate("C14", C(0,3)); - setParamNoUpdate("C15", C(0,4)); - setParamNoUpdate("C16", C(0,5)); - setParamNoUpdate("C22", C(1,1)); - setParamNoUpdate("C23", C(1,2)); - setParamNoUpdate("C24", C(1,3)); - setParamNoUpdate("C25", C(1,4)); - setParamNoUpdate("C26", C(1,5)); - setParamNoUpdate("C33", C(2,2)); - setParamNoUpdate("C34", C(2,3)); - setParamNoUpdate("C35", C(2,4)); - setParamNoUpdate("C36", C(2,5)); - setParamNoUpdate("C44", C(3,3)); - setParamNoUpdate("C45", C(3,4)); - setParamNoUpdate("C46", C(3,5)); - setParamNoUpdate("C55", C(4,4)); - setParamNoUpdate("C56", C(4,5)); - setParamNoUpdate("C66", C(5,5)); + this->Cprime = C; // set internal Cijkl matrix expressed in the canonical frame of reference this->updateInternalParameters(); - Matrix tangent(6,6); + Matrix tangent(6, 6); this->computeTangentModuliOnQuad(tangent); Real tangent_error = (tangent - C).norm(); ASSERT_NEAR(tangent_error, 0, 1e-14); } /* -------------------------------------------------------------------------- */ namespace { template class TestElasticMaterialFixture : public ::TestMaterialFixture {}; TYPED_TEST_CASE(TestElasticMaterialFixture, types); -TYPED_TEST(TestElasticMaterialFixture, ElasticComputeStress) { +TYPED_TEST(TestElasticMaterialFixture, ComputeStress) { this->material->testComputeStress(); } -TYPED_TEST(TestElasticMaterialFixture, ElasticEnergyDensity) { +TYPED_TEST(TestElasticMaterialFixture, EnergyDensity) { this->material->testEnergyDensity(); } -TYPED_TEST(TestElasticMaterialFixture, ElasticComputeTangentModuli) { +TYPED_TEST(TestElasticMaterialFixture, ComputeTangentModuli) { this->material->testComputeTangentModuli(); } -TYPED_TEST(TestElasticMaterialFixture, ElasticComputePushWaveSpeed) { +TYPED_TEST(TestElasticMaterialFixture, ComputePushWaveSpeed) { this->material->testPushWaveSpeed(); } -TYPED_TEST(TestElasticMaterialFixture, ElasticComputeShearWaveSpeed) { +TYPED_TEST(TestElasticMaterialFixture, ComputeShearWaveSpeed) { this->material->testShearWaveSpeed(); } } // namespace diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_finite_def_materials.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_finite_def_materials.cc index f1304ac25..e00b46133 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/test_finite_def_materials.cc +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_finite_def_materials.cc @@ -1,61 +1,59 @@ /* -------------------------------------------------------------------------- */ #include "material_neohookean.hh" #include "solid_mechanics_model.hh" #include "test_material_fixtures.hh" #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; using types = ::testing::Types< Traits, Traits, Traits>; /*****************************************************************/ template <> void FriendMaterial>::testComputeStress() { AKANTU_TO_IMPLEMENT(); } /*****************************************************************/ template <> void FriendMaterial>::testComputeTangentModuli() { AKANTU_TO_IMPLEMENT(); } /*****************************************************************/ template <> void FriendMaterial>::testEnergyDensity() { AKANTU_TO_IMPLEMENT(); } /*****************************************************************/ namespace { template class TestFiniteDefMaterialFixture : public ::TestMaterialFixture {}; TYPED_TEST_CASE(TestFiniteDefMaterialFixture, types); -TYPED_TEST(TestFiniteDefMaterialFixture, FiniteDefComputeStress) { +TYPED_TEST(TestFiniteDefMaterialFixture, DISABLED_ComputeStress) { this->material->testComputeStress(); } -TYPED_TEST(TestFiniteDefMaterialFixture, FiniteDefEnergyDensity) { +TYPED_TEST(TestFiniteDefMaterialFixture, DISABLED_EnergyDensity) { this->material->testEnergyDensity(); } -TYPED_TEST(TestFiniteDefMaterialFixture, FiniteDefComputeTangentModuli) { +TYPED_TEST(TestFiniteDefMaterialFixture, DISABLED_ComputeTangentModuli) { this->material->testComputeTangentModuli(); } - -TYPED_TEST(TestFiniteDefMaterialFixture, FiniteDefComputePushWaveSpeed) { +TYPED_TEST(TestFiniteDefMaterialFixture, DISABLED_ComputePushWaveSpeed) { this->material->testPushWaveSpeed(); } - -TYPED_TEST(TestFiniteDefMaterialFixture, FiniteDefComputeShearWaveSpeed) { +TYPED_TEST(TestFiniteDefMaterialFixture, DISABLED_ComputeShearWaveSpeed) { this->material->testShearWaveSpeed(); } } /*****************************************************************/ diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_fixtures.hh b/test/test_model/test_solid_mechanics_model/test_materials/test_material_fixtures.hh index f3fd93472..9eed7866b 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/test_material_fixtures.hh +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_fixtures.hh @@ -1,184 +1,170 @@ /* -------------------------------------------------------------------------- */ #include "material_elastic.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; /* -------------------------------------------------------------------------- */ template class FriendMaterial : public T { public: ~FriendMaterial() = default; virtual void testComputeStress() { AKANTU_TO_IMPLEMENT(); }; virtual void testComputeTangentModuli() { AKANTU_TO_IMPLEMENT(); }; virtual void testEnergyDensity() { AKANTU_TO_IMPLEMENT(); }; virtual void testPushWaveSpeed() { AKANTU_TO_IMPLEMENT(); } virtual void testShearWaveSpeed() { AKANTU_TO_IMPLEMENT(); } static inline Matrix getDeviatoricStrain(Real intensity); static inline Matrix getHydrostaticStrain(Real intensity); static inline const Matrix reverseRotation(Matrix mat, Matrix rotation_matrix) { Matrix R = rotation_matrix; Matrix m2 = mat * R; Matrix m1 = R.transpose(); return m1 * m2; }; static inline const Matrix applyRotation(Matrix mat, const Matrix rotation_matrix) { Matrix R = rotation_matrix; Matrix m2 = mat * R.transpose(); Matrix m1 = R; return m1 * m2; }; static inline Matrix getRandomRotation3d(); static inline Matrix getRandomRotation2d(); using T::T; }; /* -------------------------------------------------------------------------- */ template Matrix FriendMaterial::getDeviatoricStrain(Real intensity) { Matrix dev = {{0., 1., 2.}, {1., 0., 3.}, {2., 3., 0.}}; dev *= intensity; return dev; } /* -------------------------------------------------------------------------- */ template Matrix FriendMaterial::getHydrostaticStrain(Real intensity) { Matrix dev = {{1., 0., 0.}, {0., 2., 0.}, {0., 0., 3.}}; dev *= intensity; return dev; } /* -------------------------------------------------------------------------- */ template Matrix FriendMaterial::getRandomRotation3d() { constexpr UInt Dim = 3; // Seed with a real random value, if available std::random_device rd; std::mt19937 gen(rd()); // Standard mersenne_twister_engine seeded with rd() - std::uniform_real_distribution<> dis; + std::uniform_real_distribution dis; Vector v1(Dim); Vector v2(Dim); Vector v3(Dim); for (UInt i = 0; i < Dim; ++i) { v1(i) = dis(gen); v2(i) = dis(gen); } v1.normalize(); - v2.normalize(); auto d = v1.dot(v2); v2 -= v1 * d; v2.normalize(); - d = v1.dot(v2); - v2 -= v1 * d; - v2.normalize(); - v3.crossProduct(v1, v2); - Vector test_axis(3); - test_axis.crossProduct(v1, v2); - test_axis -= v3; - if (test_axis.norm() > 8 * std::numeric_limits::epsilon()) { - AKANTU_ERROR("The axis vectors do not form a right-handed coordinate " - << "system. I. e., ||n1 x n2 - n3|| should be zero, but " - << "it is " << test_axis.norm() << "."); - } - Matrix mat(Dim, Dim); for (UInt i = 0; i < Dim; ++i) { mat(0, i) = v1[i]; mat(1, i) = v2[i]; mat(2, i) = v3[i]; } return mat; } /* -------------------------------------------------------------------------- */ template Matrix FriendMaterial::getRandomRotation2d() { constexpr UInt Dim = 2; // Seed with a real random value, if available std::random_device rd; std::mt19937 gen(rd()); // Standard mersenne_twister_engine seeded with rd() std::uniform_real_distribution<> dis; // v1 and v2 are such that they form a right-hand basis with prescribed v3, // it's need (at least) for 2d linear elastic anisotropic and (orthotropic) // materials to rotate the tangent stiffness Vector v1(3); Vector v2(3); Vector v3 = {0, 0, 1}; for (UInt i = 0; i < Dim; ++i) { v1(i) = dis(gen); // v2(i) = dis(gen); } v1.normalize(); v2.crossProduct(v3, v1); Matrix mat(Dim, Dim); for (UInt i = 0; i < Dim; ++i) { mat(0, i) = v1[i]; mat(1, i) = v2[i]; } return mat; } /* -------------------------------------------------------------------------- */ template class TestMaterialBaseFixture : public ::testing::Test { public: using mat_class = typename T::mat_class; void SetUp() override { constexpr auto spatial_dimension = T::Dim; mesh = std::make_unique(spatial_dimension); model = std::make_unique(*mesh); material = std::make_unique(*model); } void TearDown() override { material.reset(nullptr); model.reset(nullptr); mesh.reset(nullptr); } using friend_class = FriendMaterial; protected: std::unique_ptr mesh; std::unique_ptr model; std::unique_ptr material; }; /* -------------------------------------------------------------------------- */ template