diff --git a/examples/implicit/implicit_dynamic.cc b/examples/implicit/implicit_dynamic.cc index 417597edf..6dd55d316 100644 --- a/examples/implicit/implicit_dynamic.cc +++ b/examples/implicit/implicit_dynamic.cc @@ -1,148 +1,148 @@ /** * @file implicit_dynamic.cc * * @author Nicolas Richart * * @date creation: Sun Oct 19 2014 * * @brief This code refers to the implicit dynamic example from the user manual * * @section LICENSE * * Copyright (©) 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 "communicator.hh" #include "non_linear_solver.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; /* -------------------------------------------------------------------------- */ const Real bar_length = 10.; const Real bar_height = 1.; const Real bar_depth = 1.; const Real F = 5e3; const Real L = bar_length; const Real I = bar_depth * bar_height * bar_height * bar_height / 12.; const Real E = 12e7; const Real rho = 1000; const Real m = rho * bar_height * bar_depth; static Real w(UInt n) { return n * n * M_PI * M_PI / (L * L) * sqrt(E * I / m); } static Real analytical_solution(Real time) { return 2 * F * L * L * L / (pow(M_PI, 4) * E * I) * ((1. - cos(w(1) * time)) + (1. - cos(w(3) * time)) / 81. + (1. - cos(w(5) * time)) / 625.); } const UInt spatial_dimension = 2; const Real time_step = 1e-4; const Real max_time = 0.62; /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize("material_dynamic.dat", argc, argv); Mesh mesh(spatial_dimension); const auto & comm = Communicator::getStaticCommunicator(); Int prank = comm.whoAmI(); if (prank == 0) mesh.read("beam.msh"); mesh.distribute(); SolidMechanicsModel model(mesh); /// model initialization model.initFull(_analysis_method = _implicit_dynamic); Material & mat = model.getMaterial(0); mat.setParam("E", E); mat.setParam("rho", rho); Array & force = model.getForce(); Array & displacment = model.getDisplacement(); // boundary conditions model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "blocked"); model.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "blocked"); model.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "roller"); const Array & trac_nodes = mesh.getElementGroup("traction").getNodes(); bool dump_node = false; if (trac_nodes.size() > 0 && mesh.isLocalOrMasterNode(trac_nodes(0))) { force(trac_nodes(0), 1) = F; dump_node = true; } // output setup std::ofstream pos; pos.open("position.csv"); if (!pos.good()) - AKANTU_DEBUG_ERROR("Cannot open file \"position.csv\""); + AKANTU_ERROR("Cannot open file \"position.csv\""); pos << "id,time,position,solution" << std::endl; model.setBaseName("dynamic"); model.addDumpFieldVector("displacement"); model.addDumpField("velocity"); model.addDumpField("acceleration"); model.addDumpField("force"); model.addDumpField("residual"); model.dump(); model.setTimeStep(time_step); auto & solver = model.getNonLinearSolver(); solver.set("max_iterations", 100); solver.set("threshold", 1e-12); solver.set("convergence_type", _scc_solution); /// time loop Real time = 0.; for (UInt s = 1; time < max_time; ++s, time += time_step) { if (prank == 0) std::cout << s << "\r" << std::flush; model.solveStep(); if (dump_node) pos << s << "," << time << "," << displacment(trac_nodes(0), 1) << "," << analytical_solution(s * time_step) << std::endl; if (s % 100 == 0) model.dump(); } std::cout << std::endl; pos.close(); finalize(); return EXIT_SUCCESS; } diff --git a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc index 27306a2ed..c8d7ac833 100644 --- a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc +++ b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc @@ -1,553 +1,553 @@ /** * @file solid_mechanics_model_RVE.cc * @author Aurelia Isabel Cuba Ramos * @date Wed Jan 13 15:32:35 2016 * * @brief Implementation of SolidMechanicsModelRVE * * @section LICENSE * * Copyright (©) 2010-2011 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_RVE.hh" #include "element_group.hh" #include "material_damage_iterative.hh" #include "node_group.hh" #include "non_linear_solver.hh" #include "parser.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ SolidMechanicsModelRVE::SolidMechanicsModelRVE(Mesh & mesh, bool use_RVE_mat_selector, UInt nb_gel_pockets, UInt dim, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, dim, id, memory_id), volume(0.), use_RVE_mat_selector(use_RVE_mat_selector), nb_gel_pockets(nb_gel_pockets), nb_dumps(0) { AKANTU_DEBUG_IN(); /// find the four corner nodes of the RVE findCornerNodes(); /// remove the corner nodes from the surface node groups: /// This most be done because corner nodes a not periodic mesh.getElementGroup("top").removeNode(corner_nodes(2)); mesh.getElementGroup("top").removeNode(corner_nodes(3)); mesh.getElementGroup("left").removeNode(corner_nodes(3)); mesh.getElementGroup("left").removeNode(corner_nodes(0)); mesh.getElementGroup("bottom").removeNode(corner_nodes(1)); mesh.getElementGroup("bottom").removeNode(corner_nodes(0)); mesh.getElementGroup("right").removeNode(corner_nodes(2)); mesh.getElementGroup("right").removeNode(corner_nodes(1)); const auto & bottom = mesh.getElementGroup("bottom").getNodeGroup(); bottom_nodes.insert(bottom.begin(), bottom.end()); const auto & left = mesh.getElementGroup("left").getNodeGroup(); left_nodes.insert(left.begin(), left.end()); // /// enforce periodicity on the displacement fluctuations // auto surface_pair_1 = std::make_pair("top", "bottom"); // auto surface_pair_2 = std::make_pair("right", "left"); // SurfacePairList surface_pairs_list; // surface_pairs_list.push_back(surface_pair_1); // surface_pairs_list.push_back(surface_pair_2); // TODO: To Nicolas correct the PBCs // this->setPBC(surface_pairs_list); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelRVE::~SolidMechanicsModelRVE() = default; /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); auto options_cp(options); options_cp.analysis_method = AnalysisMethod::_static; SolidMechanicsModel::initFullImpl(options_cp); this->initMaterials(); auto & fem = this->getFEEngine("SolidMechanicsFEEngine"); /// compute the volume of the RVE GhostType gt = _not_ghost; for (auto element_type : this->mesh.elementTypes(spatial_dimension, gt, _ek_not_defined)) { Array Volume(this->mesh.getNbElement(element_type) * fem.getNbIntegrationPoints(element_type), 1, 1.); this->volume = fem.integrate(Volume, element_type); } std::cout << "The volume of the RVE is " << this->volume << std::endl; /// dumping std::stringstream base_name; base_name << this->id; // << this->memory_id - 1; this->setBaseName(base_name.str()); this->addDumpFieldVector("displacement"); this->addDumpField("stress"); this->addDumpField("grad_u"); this->addDumpField("eigen_grad_u"); this->addDumpField("blocked_dofs"); this->addDumpField("material_index"); this->addDumpField("damage"); this->addDumpField("Sc"); this->addDumpField("external_force"); this->addDumpField("equivalent_stress"); this->addDumpField("internal_force"); this->dump(); this->nb_dumps += 1; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::applyBoundaryConditions( const Matrix & displacement_gradient) { AKANTU_DEBUG_IN(); /// get the position of the nodes const Array & pos = mesh.getNodes(); /// storage for the coordinates of a given node and the displacement that will /// be applied Vector x(spatial_dimension); Vector appl_disp(spatial_dimension); /// fix top right node UInt node = this->corner_nodes(2); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); // (*this->blocked_dofs)(node,0) = true; (*this->displacement)(node,0) = 0.; // (*this->blocked_dofs)(node,1) = true; (*this->displacement)(node,1) = 0.; /// apply Hx at all the other corner nodes; H: displ. gradient node = this->corner_nodes(0); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); node = this->corner_nodes(1); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); node = this->corner_nodes(3); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::findCornerNodes() { AKANTU_DEBUG_IN(); // find corner nodes const auto & position = mesh.getNodes(); const auto & lower_bounds = mesh.getLowerBounds(); const auto & upper_bounds = mesh.getUpperBounds(); AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "This is 2D only!"); corner_nodes.resize(4); corner_nodes.set(UInt(-1)); for (auto && data : enumerate(make_view(position, spatial_dimension))) { auto node = std::get<0>(data); const auto & X = std::get<1>(data); auto distance = X.distance(lower_bounds); // node 1 if (Math::are_float_equal(distance, 0)) { corner_nodes(0) = node; } // node 2 else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && Math::are_float_equal(X(_y), lower_bounds(_y))) { corner_nodes(1) = node; } // node 3 else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && Math::are_float_equal(X(_y), upper_bounds(_y))) { corner_nodes(2) = node; } // node 4 else if (Math::are_float_equal(X(_x), lower_bounds(_x)) && Math::are_float_equal(X(_y), upper_bounds(_y))) { corner_nodes(3) = node; } } for (UInt i = 0; i < corner_nodes.size(); ++i) { if (corner_nodes(i) == UInt(-1)) - AKANTU_DEBUG_ERROR("The corner node " << i + 1 << " wasn't found"); + AKANTU_ERROR("The corner node " << i + 1 << " wasn't found"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::advanceASR(const Matrix & prestrain) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "This is 2D only!"); /// apply the new eigenstrain for (auto element_type : mesh.elementTypes(_element_kind = _ek_not_defined)) { Array & prestrain_vect = const_cast &>(this->getMaterial("gel").getInternal( "eigen_grad_u")(element_type)); auto prestrain_it = prestrain_vect.begin(spatial_dimension, spatial_dimension); auto prestrain_end = prestrain_vect.end(spatial_dimension, spatial_dimension); for (; prestrain_it != prestrain_end; ++prestrain_it) (*prestrain_it) = prestrain; } /// advance the damage MaterialDamageIterative<2> & mat_paste = dynamic_cast &>(*this->materials[1]); MaterialDamageIterative<2> & mat_aggregate = dynamic_cast &>(*this->materials[0]); UInt nb_damaged_elements = 0; Real max_eq_stress_aggregate = 0; Real max_eq_stress_paste = 0; auto & solver = this->getNonLinearSolver(); solver.set("max_iterations", 2); solver.set("threshold", 1e-6); solver.set("convergence_type", _scc_solution); do { this->solveStep(); /// compute damage max_eq_stress_aggregate = mat_aggregate.getNormMaxEquivalentStress(); max_eq_stress_paste = mat_paste.getNormMaxEquivalentStress(); nb_damaged_elements = 0; if (max_eq_stress_aggregate > max_eq_stress_paste) nb_damaged_elements = mat_aggregate.updateDamage(); else if (max_eq_stress_aggregate < max_eq_stress_paste) nb_damaged_elements = mat_paste.updateDamage(); else nb_damaged_elements = (mat_paste.updateDamage() + mat_aggregate.updateDamage()); std::cout << "the number of damaged elements is " << nb_damaged_elements << std::endl; } while (nb_damaged_elements); if (this->nb_dumps % 10 == 0) { this->dump(); } this->nb_dumps += 1; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModelRVE::averageTensorField(UInt row_index, UInt col_index, const ID & field_type) { AKANTU_DEBUG_IN(); auto & fem = this->getFEEngine("SolidMechanicsFEEngine"); Real average = 0; for (auto element_type : mesh.elementTypes(_element_kind = _ek_not_defined)) { if (field_type == "stress") { for (UInt m = 0; m < this->materials.size(); ++m) { const auto & stress_vec = this->materials[m]->getStress(element_type); const auto & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_stress_vec(elem_filter.size(), spatial_dimension * spatial_dimension, "int_of_stress"); fem.integrate(stress_vec, int_stress_vec, spatial_dimension * spatial_dimension, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) average += int_stress_vec( k, row_index * spatial_dimension + col_index); // 3 is the value // for the yy (in // 3D, the value is // 4) } } else if (field_type == "strain") { for (UInt m = 0; m < this->materials.size(); ++m) { const auto & gradu_vec = this->materials[m]->getGradU(element_type); const auto & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_gradu_vec(elem_filter.size(), spatial_dimension * spatial_dimension, "int_of_gradu"); fem.integrate(gradu_vec, int_gradu_vec, spatial_dimension * spatial_dimension, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) /// averaging is done only for normal components, so stress and strain /// are equal average += 0.5 * (int_gradu_vec(k, row_index * spatial_dimension + col_index) + int_gradu_vec(k, col_index * spatial_dimension + row_index)); } } else if (field_type == "eigen_grad_u") { for (UInt m = 0; m < this->materials.size(); ++m) { const auto & eigen_gradu_vec = this->materials[m]->getInternal("eigen_grad_u")(element_type); const auto & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_eigen_gradu_vec(elem_filter.size(), spatial_dimension * spatial_dimension, "int_of_gradu"); fem.integrate(eigen_gradu_vec, int_eigen_gradu_vec, spatial_dimension * spatial_dimension, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) /// averaging is done only for normal components, so stress and strain /// are equal average += int_eigen_gradu_vec(k, row_index * spatial_dimension + col_index); } } else { - AKANTU_DEBUG_ERROR("Averaging not implemented for this field!!!"); + AKANTU_ERROR("Averaging not implemented for this field!!!"); } } return average / this->volume; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::homogenizeStiffness(Matrix & C_macro) { AKANTU_DEBUG_IN(); const UInt dim = 2; AKANTU_DEBUG_ASSERT(this->spatial_dimension == dim, "Is only implemented for 2D!!!"); /// apply three independent loading states to determine C /// 1. eps_el = (1;0;0) 2. eps_el = (0,1,0) 3. eps_el = (0,0,0.5) /// clear the eigenstrain Matrix zero_eigengradu(dim, dim, 0.); GhostType gt = _not_ghost; for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { auto & prestrain_vect = const_cast &>(this->getMaterial("gel").getInternal( "eigen_grad_u")(element_type)); auto prestrain_it = prestrain_vect.begin(spatial_dimension, spatial_dimension); auto prestrain_end = prestrain_vect.end(spatial_dimension, spatial_dimension); for (; prestrain_it != prestrain_end; ++prestrain_it) (*prestrain_it) = zero_eigengradu; } /// storage for results of 3 different loading states UInt voigt_size = VoigtHelper::size; Matrix stresses(voigt_size, voigt_size, 0.); Matrix strains(voigt_size, voigt_size, 0.); Matrix H(dim, dim, 0.); /// save the damage state before filling up cracks // ElementTypeMapReal saved_damage("saved_damage"); // saved_damage.initialize(getFEEngine(), _nb_component = 1, _default_value = // 0); // this->fillCracks(saved_damage); /// virtual test 1: H(0, 0) = 0.01; this->performVirtualTesting(H, stresses, strains, 0); /// virtual test 2: H.clear(); H(1, 1) = 0.01; this->performVirtualTesting(H, stresses, strains, 1); /// virtual test 3: H.clear(); H(0, 1) = 0.01; this->performVirtualTesting(H, stresses, strains, 2); /// drain cracks // this->drainCracks(saved_damage); /// compute effective stiffness Matrix eps_inverse(voigt_size, voigt_size); eps_inverse.inverse(strains); C_macro.mul(stresses, eps_inverse); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::performVirtualTesting(const Matrix & H, Matrix & eff_stresses, Matrix & eff_strains, const UInt test_no) { AKANTU_DEBUG_IN(); this->applyBoundaryConditions(H); auto & solver = this->getNonLinearSolver(); solver.set("max_iterations", 2); solver.set("threshold", 1e-6); solver.set("convergence_type", _scc_solution); this->solveStep(); /// get average stress and strain eff_stresses(0, test_no) = this->averageTensorField(0, 0, "stress"); eff_strains(0, test_no) = this->averageTensorField(0, 0, "strain"); eff_stresses(1, test_no) = this->averageTensorField(1, 1, "stress"); eff_strains(1, test_no) = this->averageTensorField(1, 1, "strain"); eff_stresses(2, test_no) = this->averageTensorField(1, 0, "stress"); eff_strains(2, test_no) = 2. * this->averageTensorField(1, 0, "strain"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::homogenizeEigenGradU( Matrix & eigen_gradu_macro) { AKANTU_DEBUG_IN(); eigen_gradu_macro(0, 0) = this->averageTensorField(0, 0, "eigen_grad_u"); eigen_gradu_macro(1, 1) = this->averageTensorField(1, 1, "eigen_grad_u"); eigen_gradu_macro(0, 1) = this->averageTensorField(0, 1, "eigen_grad_u"); eigen_gradu_macro(1, 0) = this->averageTensorField(1, 0, "eigen_grad_u"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initMaterials() { AKANTU_DEBUG_IN(); // make sure the material are instantiated if (!are_materials_instantiated) instantiateMaterials(); if (use_RVE_mat_selector) { const Vector & lowerBounds = mesh.getLowerBounds(); const Vector & upperBounds = mesh.getUpperBounds(); Real bottom = lowerBounds(1); Real top = upperBounds(1); Real box_size = std::abs(top - bottom); Real eps = box_size * 1e-6; auto tmp = std::make_shared(*this, box_size, "gel", this->nb_gel_pockets, eps); tmp->setFallback(material_selector); material_selector = tmp; } SolidMechanicsModel::initMaterials(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::fillCracks(ElementTypeMapReal & saved_damage) { const auto & mat_gel = this->getMaterial("gel"); Real E_gel = mat_gel.get("E"); Real E_homogenized = 0.; for (auto && mat : materials) { if (mat->getName() == "gel" || mat->getName() == "FE2_mat") continue; Real E = mat->get("E"); auto & damage = mat->getInternal("damage"); for (auto && type : damage.elementTypes()) { const auto & elem_filter = mat->getElementFilter(type); auto nb_integration_point = getFEEngine().getNbIntegrationPoints(type); auto sav_dam_it = make_view(saved_damage(type), nb_integration_point).begin(); for (auto && data : zip(elem_filter, make_view(damage(type), nb_integration_point))) { auto el = std::get<0>(data); auto & dam = std::get<1>(data); Vector sav_dam = sav_dam_it[el]; sav_dam = dam; for (auto q : arange(dam.size())) { E_homogenized = (E_gel - E) * dam(q) + E; dam(q) = 1. - (E_homogenized / E); } } } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::drainCracks( const ElementTypeMapReal & saved_damage) { for (auto && mat : materials) { if (mat->getName() == "gel" || mat->getName() == "FE2_mat") continue; auto & damage = mat->getInternal("damage"); for (auto && type : damage.elementTypes()) { const auto & elem_filter = mat->getElementFilter(type); auto nb_integration_point = getFEEngine().getNbIntegrationPoints(type); auto sav_dam_it = make_view(saved_damage(type), nb_integration_point).begin(); for (auto && data : zip(elem_filter, make_view(damage(type), nb_integration_point))) { auto el = std::get<0>(data); auto & dam = std::get<1>(data); Vector sav_dam = sav_dam_it[el]; dam = sav_dam; } } } } } // namespace akantu diff --git a/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_inline_impl.cc b/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_inline_impl.cc index 201eb94bc..8b15de85b 100644 --- a/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_inline_impl.cc +++ b/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_inline_impl.cc @@ -1,177 +1,177 @@ /** * @file material_vreepeerlings_inline_impl.cc * * @author Cyprien Wolff * * * @brief Specialization of the material class for the VreePeerlings material * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template class MatParent> inline void MaterialVreePeerlings::computeStressOnQuad(Matrix & grad_u, Matrix & sigma, Real & dam, Real & Equistrain, Real & Equistrain_rate, Real & Kapaq, Real dt, Matrix & strain_rate_vrplgs, Real & FullDam_ValStrain, Real & FullDam_ValStrain_rate, Real & Nb_damage) { Real I1=0.; /// trace initialization of the strain tensor Real J2=0.; /// J2 [J2=1/6(3 strain:strain - I1*I1)] initialization Real I1rate=0.; /// trace initialization of the strain rate tensor Real J2rate=0.; /// J2 [J2=1/3(3 strain:strainrate - I1*I1rate)] initialization if (spatial_dimension == 1){ I1 = grad_u.trace(); J2 = .5*grad_u(0,0)*grad_u(0,0) - I1 * I1/6.; I1rate = strain_rate_vrplgs.trace(); J2rate = grad_u(0,0)*strain_rate_vrplgs(0,0) -I1 * I1rate/6.; } else { // if(this->plane_stress) { // Real tmp = this->nu/(this->nu - 1); // tmp *= tmp; // I1 = (grad_u(0,0) + grad_u(1,1))*(1 - 2*this->nu)/(1 - this->nu); // J2 = .5*(grad_u(0,0)*grad_u(0,0) + // grad_u(1,1)*grad_u(1,1) + // tmp*(grad_u(0,0) + grad_u(1,1))*(grad_u(0,0) + grad_u(1,1)) + // .5*(grad_u(0,1) + grad_u(1,0))*(grad_u(0,1) + grad_u(1,0))) - // I1*I1/6.; // I1rate = (strain_rate_vrplgs(0,0) + strain_rate_vrplgs(1,1))*(1 - 2*this->nu)/(1 - this->nu); // J2rate = (grad_u(0,0)*strain_rate_vrplgs(0,0) + // grad_u(1,1)*strain_rate_vrplgs(1,1) + // tmp*(grad_u(0,0) + grad_u(1,1))*(strain_rate_vrplgs(0,0) + strain_rate_vrplgs(1,1)) + // (grad_u(0,1)*strain_rate_vrplgs(1,0)) + (grad_u(0,1)*strain_rate_vrplgs(1,0))) - // I1*I1rate/3.; // } // else { I1 = grad_u.trace(); for(UInt i=0; inu); Real tmp_2 = (12 * Kct)/((1 + this->nu)*(1 + this->nu)); Equistrain = tmp_1 * I1 / (2*Kct) + 1./(2*Kct) * std::sqrt(tmp_1*tmp_1*I1*I1 + tmp_2*J2); if(I1 * I1rate > 0 || J2rate > 0){ Equistrain_rate =tmp_1 * std::abs(I1rate) / (2*Kct) + 1./(4*Kct)*(2*tmp_1*tmp_1*I1*I1rate + tmp_2*J2rate) / std::sqrt(tmp_1*tmp_1*I1*I1 + tmp_2*J2); }else{ - AKANTU_DEBUG_ERROR("This instruction here was commented but has to be checked"); + AKANTU_ERROR("This instruction here was commented but has to be checked"); //Equistrain_rate = Equistrain_rate; } if(!this->is_non_local) { computeDamageAndStressOnQuad(sigma, dam, Equistrain, Equistrain_rate, Kapaq, dt, FullDam_ValStrain, FullDam_ValStrain_rate, Nb_damage); } } /* -------------------------------------------------------------------------- */ template class MatParent> inline void MaterialVreePeerlings::computeDamageAndStressOnQuad(Matrix & sigma, Real & dam, Real & Equistrain, Real & Equistrain_rate, Real & Kapaq, Real dt, Real & FullDam_ValStrain, Real & FullDam_ValStrain_rate, Real & Nb_damage) { //--------------------------------------------------------------------------------------- // rate-dependence model // //--------------------------------------------------------------------------------------- Real Kapao = Crate*(1.+(std::pow(std::abs(Equistrain_rate)*Arate,Brate)))*(1.-Kapao_randomness) + Kapao_randomness*Kapaq; Real Kapac_99 = Kapac*.99; Real Kapaodyn = 0; if (Kapao > Kapaoi) { if (Kapao < Kapac_99) { Kapaodyn = Kapao; } else { Kapaodyn = Kapac_99; } } else { Kapaodyn = Kapaoi; } Real Fd1 = Equistrain - Kapaoi; if (Fd1 > 0) { Real dam1 = 1. - Kapaodyn/Equistrain * ((Kapac - Equistrain)/(Kapac - Kapaodyn)); if (dam1 > dam){ dam = std::min(dam1,1.); Nb_damage = Nb_damage + 1.; if (dam >= 1.){ FullDam_ValStrain = Equistrain; FullDam_ValStrain_rate = Equistrain_rate; } } } //--------------------------------------------------------------------------------------- // delayed damage (see Marions thesis page 68) //--------------------------------------------------------------------------------------- // Real viscosity = 10.; // Real damRateInfini = 10000000.; // Real gequi = 1. - Kapaq/Equistrain * ((Kapac-Equistrain)/(Kapac - Kapaq)); // if (gequi - dam > 0){ // // Real damRate = damRateInfini * (1. - std::exp(-1*viscosity * (gequi - dam))); // if (damrate > 0){ // // if (dam < 1.){ // dam = dam + damRate*dt; // } else { // dam = 1.; // } // } // } // //--------------------------------------------------------------------------------------- sigma *= 1-dam; } /* -------------------------------------------------------------------------- */ diff --git a/extra_packages/parallel-cohesive-element/test/test_data_accessor.hh b/extra_packages/parallel-cohesive-element/test/test_data_accessor.hh index ec366590e..7ac9bdd60 100644 --- a/extra_packages/parallel-cohesive-element/test/test_data_accessor.hh +++ b/extra_packages/parallel-cohesive-element/test/test_data_accessor.hh @@ -1,134 +1,134 @@ /** * @file test_data_accessor.hh * * @author Marco Vocialta * @author Nicolas Richart * * @date creation: Thu Apr 11 2013 * @date last modification: Thu Jun 05 2014 * * @brief Data Accessor class for testing * * @section LICENSE * * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "data_accessor.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class TestAccessor : public DataAccessor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: inline TestAccessor(const Mesh & mesh, const ElementTypeMapArray & barycenters); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Barycenter, barycenters, Real); /* ------------------------------------------------------------------------ */ /* Ghost Synchronizer inherited members */ /* ------------------------------------------------------------------------ */ protected: inline UInt getNbDataForElements(const Array & elements, SynchronizationTag tag) const; inline void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const; inline void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: const ElementTypeMapArray & barycenters; const Mesh & mesh; }; /* -------------------------------------------------------------------------- */ /* TestSynchronizer implementation */ /* -------------------------------------------------------------------------- */ inline TestAccessor::TestAccessor(const Mesh & mesh, const ElementTypeMapArray & barycenters) : barycenters(barycenters), mesh(mesh) { } inline UInt TestAccessor::getNbDataForElements(const Array & elements, __attribute__ ((unused)) SynchronizationTag tag) const { if(elements.getSize()) // return Mesh::getSpatialDimension(elements(0).type) * sizeof(Real) * elements.getSize(); return mesh.getSpatialDimension() * sizeof(Real) * elements.getSize(); else return 0; } inline void TestAccessor::packElementData(CommunicationBuffer & buffer, const Array & elements, __attribute__ ((unused)) SynchronizationTag tag) const { UInt spatial_dimension = mesh.getSpatialDimension(); Array::const_iterator bit = elements.begin(); Array::const_iterator bend = elements.end(); for (; bit != bend; ++bit) { const Element & element = *bit; Vector bary(this->barycenters(element.type, element.ghost_type).storage() + element.element * spatial_dimension, spatial_dimension); buffer << bary; } } inline void TestAccessor::unpackElementData(CommunicationBuffer & buffer, const Array & elements, __attribute__ ((unused)) SynchronizationTag tag) { UInt spatial_dimension = mesh.getSpatialDimension(); Array::const_iterator bit = elements.begin(); Array::const_iterator bend = elements.end(); for (; bit != bend; ++bit) { const Element & element = *bit; Vector barycenter_loc(this->barycenters(element.type, element.ghost_type).storage() + element.element * spatial_dimension, spatial_dimension); Vector bary(spatial_dimension); buffer >> bary; std::cout << element << barycenter_loc << std::endl; Real tolerance = 1e-15; for (UInt i = 0; i < spatial_dimension; ++i) { if(!(std::abs(bary(i) - barycenter_loc(i)) <= tolerance)) - AKANTU_DEBUG_ERROR("Unpacking an unknown value for the element: " + AKANTU_ERROR("Unpacking an unknown value for the element: " << element << "(barycenter = " << barycenter_loc << " and buffer = " << bary << ") direction (" << i << ")- tag: " << tag); } } } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/functions/boundary_functions.cc b/extra_packages/traction-at-split-node-contact/src/functions/boundary_functions.cc index 2b2185404..0594b0743 100644 --- a/extra_packages/traction-at-split-node-contact/src/functions/boundary_functions.cc +++ b/extra_packages/traction-at-split-node-contact/src/functions/boundary_functions.cc @@ -1,63 +1,63 @@ /** * @file boundary_functions.cc * * * * @brief * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ #include "boundary_functions.hh" namespace akantu { /* -------------------------------------------------------------------------- */ Real integrateResidual(const std::string & sub_boundary_name, const SolidMechanicsModel & model, UInt dir) { Real int_res = 0.; const Mesh & mesh = model.getMesh(); const Array & residual = model.getResidual(); // do not need try catch, as all subboundaries should be everywhere. // try { const ElementGroup & boundary = mesh.getElementGroup(sub_boundary_name); ElementGroup::const_node_iterator nit = boundary.node_begin(); ElementGroup::const_node_iterator nend = boundary.node_end(); for (; nit != nend; ++nit) { bool is_local_node = mesh.isLocalOrMasterNode(*nit); if (is_local_node) { int_res += residual(*nit, dir); } } // } catch(debug::Exception e) { - // // AKANTU_DEBUG_ERROR("Error computing integrateResidual. Cannot get SubBoundary: " + // // AKANTU_ERROR("Error computing integrateResidual. Cannot get SubBoundary: " // // << sub_boundary_name << " [" << e.what() << "]"); // } StaticCommunicator::getStaticCommunicator().allReduce(&int_res, 1, _so_sum); return int_res; } /* -------------------------------------------------------------------------- */ void boundaryFix(Mesh & mesh, const std::vector & sub_boundary_names) { std::vector::const_iterator it = sub_boundary_names.begin(); std::vector::const_iterator end = sub_boundary_names.end(); for (; it != end; ++it) { if (mesh.element_group_find(*it) == mesh.element_group_end()) { mesh.createElementGroup(*it,mesh.getSpatialDimension()-1); // empty element group } } } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/functions/node_filter.hh b/extra_packages/traction-at-split-node-contact/src/functions/node_filter.hh index 175d71d1a..daad5b664 100644 --- a/extra_packages/traction-at-split-node-contact/src/functions/node_filter.hh +++ b/extra_packages/traction-at-split-node-contact/src/functions/node_filter.hh @@ -1,106 +1,98 @@ /** * @file node_filter.hh * * @author David Simon Kammer * * * @brief to filter nodes with functors * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #ifndef __AST_NODE_FILTER_HH__ #define __AST_NODE_FILTER_HH__ /* -------------------------------------------------------------------------- */ // akantu #include "aka_common.hh" #include "mesh_filter.hh" - namespace akantu { /* -------------------------------------------------------------------------- */ class GeometryFilter : public NodeFilterFunctor { public: - GeometryFilter(const Mesh & mesh, - UInt dir, - Real limit) : - NodeFilterFunctor(), - mesh(mesh), dir(dir), limit(limit) { + GeometryFilter(const Mesh & mesh, UInt dir, Real limit) + : NodeFilterFunctor(), mesh(mesh), dir(dir), limit(limit) { this->positions = &(mesh.getNodes()); }; - ~GeometryFilter() {}; - - bool operator() (UInt node) { - AKANTU_DEBUG_TO_IMPLEMENT(); - }; + ~GeometryFilter(){}; + + bool operator()(UInt node) { AKANTU_TO_IMPLEMENT(); }; protected: const Mesh & mesh; UInt dir; Real limit; const Array * positions; }; /* -------------------------------------------------------------------------- */ class FilterPositionsGreaterThan : public GeometryFilter { public: - FilterPositionsGreaterThan(const Mesh & mesh, - UInt dir, - Real limit) : GeometryFilter(mesh, dir, limit) {}; - ~FilterPositionsGreaterThan() {}; + FilterPositionsGreaterThan(const Mesh & mesh, UInt dir, Real limit) + : GeometryFilter(mesh, dir, limit){}; + ~FilterPositionsGreaterThan(){}; - bool operator() (UInt node) { + bool operator()(UInt node) { AKANTU_DEBUG_IN(); bool to_filter = true; - if((*this->positions)(node,this->dir) > this->limit) + if ((*this->positions)(node, this->dir) > this->limit) to_filter = false; AKANTU_DEBUG_OUT(); return to_filter; }; }; /* -------------------------------------------------------------------------- */ class FilterPositionsLessThan : public GeometryFilter { public: - FilterPositionsLessThan(const Mesh & mesh, - UInt dir, - Real limit) : GeometryFilter(mesh, dir, limit) {}; - ~FilterPositionsLessThan() {}; + FilterPositionsLessThan(const Mesh & mesh, UInt dir, Real limit) + : GeometryFilter(mesh, dir, limit){}; + ~FilterPositionsLessThan(){}; - bool operator() (UInt node) { + bool operator()(UInt node) { AKANTU_DEBUG_IN(); bool to_filter = true; - if((*this->positions)(node,this->dir) < this->limit) + if ((*this->positions)(node, this->dir) < this->limit) to_filter = false; AKANTU_DEBUG_OUT(); return to_filter; }; }; /* -------------------------------------------------------------------------- */ // this filter is erase because the convention of filter has changed!! // filter == true -> keep node // template // void applyNodeFilter(Array & nodes, FilterType & filter) { // Array::iterator<> it = nodes.begin(); // for (; it != nodes.end(); ++it) { // if (filter(*it)) { // it = nodes.erase(it); // } // } -// }; +// }; } // namespace akantu #endif /* __AST_NODE_FILTER_HH__ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.hh index a14228f58..6452eda7a 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.hh @@ -1,161 +1,161 @@ /** * @file ntn_base_friction.hh * * @author David Simon Kammer * * * @brief base class for ntn and ntrf friction * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #ifndef __AST_NTN_BASE_FRICTION_HH__ #define __AST_NTN_BASE_FRICTION_HH__ /* -------------------------------------------------------------------------- */ // akantu #include "parsable.hh" // simtools #include "ntn_base_contact.hh" /* -------------------------------------------------------------------------- */ namespace akantu { template<> inline void ParsableParamTyped< akantu::SynchronizedArray >::parseParam(const ParserParameter & in_param) { ParsableParam::parseParam(in_param); Real tmp = in_param; param.setAndChangeDefault(tmp); } /* -------------------------------------------------------------------------- */ template<> template<> inline void ParsableParamTyped< akantu::SynchronizedArray >::setTyped(const Real & value) { param.setAndChangeDefault(value); } /* -------------------------------------------------------------------------- */ class NTNBaseFriction : protected Memory, public Parsable, public Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NTNBaseFriction(NTNBaseContact * contact, const FrictionID & id = "friction", const MemoryID & memory_id = 0); virtual ~NTNBaseFriction() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// compute friction traction virtual void computeFrictionTraction(); /// compute stick traction (friction traction needed to stick the nodes) virtual void computeStickTraction(); /// apply the friction force virtual void applyFrictionTraction(); /// compute slip virtual void updateSlip(); /// register Syncronizedarrays for sync virtual void registerSynchronizedArray(SynchronizedArrayBase & array); /// dump restart file virtual void dumpRestart(const std::string & file_name) const; /// read restart file virtual void readRestart(const std::string & file_name); /// set to steady state virtual void setToSteadyState() { - AKANTU_DEBUG_TO_IMPLEMENT(); + AKANTU_TO_IMPLEMENT(); }; /// get the number of sticking nodes (in parallel) /// a node that is not in contact does not count as sticking virtual UInt getNbStickingNodes() const; /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// compute frictional strength according to friction law virtual void computeFrictionalStrength() = 0; /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Contact, contact, const NTNBaseContact *) AKANTU_GET_MACRO(IsSticking, is_sticking, const SynchronizedArray &) AKANTU_GET_MACRO(FrictionalStrength, frictional_strength, const SynchronizedArray &) AKANTU_GET_MACRO(FrictionTraction, friction_traction, const SynchronizedArray &) AKANTU_GET_MACRO(Slip, slip, const SynchronizedArray &) AKANTU_GET_MACRO(CumulativeSlip, cumulative_slip, const SynchronizedArray &) AKANTU_GET_MACRO(SlipVelocity, slip_velocity, const SynchronizedArray &) /// set parameter of a given node /// (if you need to set to all: used the setMixed function of the Parsable). virtual void setParam(const std::string & name, UInt node, Real value); // replaced by the setMixed of the Parsable // virtual void setParam(const std::string & param, Real value) { - // AKANTU_DEBUG_ERROR("Friction does not know the following parameter: " << param); + // AKANTU_ERROR("Friction does not know the following parameter: " << param); // }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: NTNBaseContact * contact; // if node is sticking SynchronizedArray is_sticking; // frictional strength SynchronizedArray frictional_strength; // friction force SynchronizedArray friction_traction; // slip SynchronizedArray slip; SynchronizedArray cumulative_slip; // slip velocity (tangential vector) SynchronizedArray slip_velocity; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "ntn_base_friction_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const NTNBaseFriction & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* __AST_NTN_BASE_FRICTION_HH__ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_initiation_function.cc b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_initiation_function.cc index 4d3932e12..ca6701f64 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_initiation_function.cc +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_initiation_function.cc @@ -1,398 +1,398 @@ /** * @file ntn_initiation_function.cc * * @author David Simon Kammer * * * @brief implementation of initializing ntn and ntrf friction * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ // simtools #include "mIIasym_contact.hh" #include "ntn_initiation_function.hh" #include "ntn_friction.hh" #include "ntrf_friction.hh" // friction regularisations #include "ntn_fricreg_rubin_ampuero.hh" #include "ntn_fricreg_simplified_prakash_clifton.hh" // friction laws #include "ntn_friclaw_linear_cohesive.hh" #include "ntn_friclaw_linear_slip_weakening.hh" #include "ntn_friclaw_linear_slip_weakening_no_healing.hh" namespace akantu { NTNBaseFriction * initializeNTNFriction(NTNBaseContact * contact, ParameterReader & data) { AKANTU_DEBUG_IN(); const std::string & friction_law = data.get("friction_law"); const std::string & friction_reg = data.get("friction_regularisation"); NTNBaseFriction * friction; bool is_ntn_contact = true; if (dynamic_cast(contact) != NULL || dynamic_cast(contact) != NULL) { is_ntn_contact = false; } if (friction_law == "coulomb") { if (friction_reg == "no_regularisation") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); } else if (friction_reg == "rubin_ampuero") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); friction->setMixed< SynchronizedArray >("t_star", data.get("t_star")); } else if (friction_reg == "simplified_prakash_clifton") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); friction->setMixed< SynchronizedArray >("t_star", data.get("t_star")); } else { - AKANTU_DEBUG_ERROR("Do not know the following friction regularisation: " + AKANTU_ERROR("Do not know the following friction regularisation: " << friction_reg); } friction->setMixed< SynchronizedArray >("mu_s", data.get("mu_s")); } // Friction Law: Linear Slip Weakening else if (friction_law == "linear_slip_weakening") { if (friction_reg == "no_regularisation") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); } else if (friction_reg == "rubin_ampuero") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); friction->setMixed< SynchronizedArray >("t_star", data.get("t_star")); } else if (friction_reg == "simplified_prakash_clifton") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); friction->setMixed< SynchronizedArray >("t_star", data.get("t_star")); } else { - AKANTU_DEBUG_ERROR("Do not know the following friction regularisation: " + AKANTU_ERROR("Do not know the following friction regularisation: " << friction_reg); } friction->setMixed< SynchronizedArray >("mu_s", data.get("mu_s")); friction->setMixed< SynchronizedArray >("mu_k", data.get("mu_k")); friction->setMixed< SynchronizedArray >("d_c", data.get("d_c")); } // Friction Law: Linear Slip Weakening No Healing else if (friction_law == "linear_slip_weakening_no_healing") { if (friction_reg == "no_regularisation") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); } else if (friction_reg == "rubin_ampuero") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); friction->setMixed< SynchronizedArray >("t_star", data.get("t_star")); } else if (friction_reg == "simplified_prakash_clifton") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); friction->setMixed< SynchronizedArray >("t_star", data.get("t_star")); } else { - AKANTU_DEBUG_ERROR("Do not know the following friction regularisation: " + AKANTU_ERROR("Do not know the following friction regularisation: " << friction_reg); } friction->setMixed< SynchronizedArray >("mu_s", data.get("mu_s")); friction->setMixed< SynchronizedArray >("mu_k", data.get("mu_k")); friction->setMixed< SynchronizedArray >("d_c", data.get("d_c")); } // Friction Law: Linear Cohesive else if (friction_law == "linear_cohesive") { if (friction_reg == "no_regularisation") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); } else if (friction_reg == "rubin_ampuero") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); friction->setMixed< SynchronizedArray >("t_star", data.get("t_star")); } else if (friction_reg == "simplified_prakash_clifton") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); friction->setMixed< SynchronizedArray >("t_star", data.get("t_star")); } else { - AKANTU_DEBUG_ERROR("Do not know the following friction regularisation: " + AKANTU_ERROR("Do not know the following friction regularisation: " << friction_reg); } friction->setMixed< SynchronizedArray >("G_c", data.get("G_c")); friction->setMixed< SynchronizedArray >("tau_c", data.get("tau_c")); friction->setMixed< SynchronizedArray >("tau_r", data.get("tau_r")); } else { - AKANTU_DEBUG_ERROR("Do not know the following friction law: " + AKANTU_ERROR("Do not know the following friction law: " << friction_law); } AKANTU_DEBUG_OUT(); return friction; } /* -------------------------------------------------------------------------- */ NTNBaseFriction * initializeNTNFriction(NTNBaseContact * contact) { AKANTU_DEBUG_IN(); std::pair sub_sect = getStaticParser().getSubSections(_st_friction); Parser::const_section_iterator it = sub_sect.first; const ParserSection & section = *it; std::string friction_law = section.getName(); std::string friction_reg = section.getOption(); if (friction_reg == "") { std::string friction_reg = "no_regularisation"; } NTNBaseFriction * friction = initializeNTNFriction(contact, friction_law, friction_reg); friction->parseSection(section); if (++it != sub_sect.second) { AKANTU_DEBUG_WARNING("There were several friction sections in input file. " << "Only first one was used and all others ignored."); } AKANTU_DEBUG_OUT(); return friction; } /* -------------------------------------------------------------------------- */ NTNBaseFriction * initializeNTNFriction(NTNBaseContact * contact, const std::string & friction_law, const std::string & friction_reg) { AKANTU_DEBUG_IN(); NTNBaseFriction * friction; // check whether is is node-to-rigid-flat contact or mIIasym (which is also ntrf) bool is_ntn_contact = true; if (dynamic_cast(contact) != NULL || dynamic_cast(contact) != NULL) { is_ntn_contact = false; } if (friction_law == "coulomb") { if (friction_reg == "no_regularisation") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); } else if (friction_reg == "rubin_ampuero") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); } else if (friction_reg == "simplified_prakash_clifton") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); } else { - AKANTU_DEBUG_ERROR("Do not know the following friction regularisation: " + AKANTU_ERROR("Do not know the following friction regularisation: " << friction_reg); } } // Friction Law: Linear Slip Weakening else if (friction_law == "linear_slip_weakening") { if (friction_reg == "no_regularisation") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); } else if (friction_reg == "rubin_ampuero") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); } else if (friction_reg == "simplified_prakash_clifton") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); } else { - AKANTU_DEBUG_ERROR("Do not know the following friction regularisation: " + AKANTU_ERROR("Do not know the following friction regularisation: " << friction_reg); } } // Friction Law: Linear Slip Weakening No Healing else if (friction_law == "linear_slip_weakening_no_healing") { if (friction_reg == "no_regularisation") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); } else if (friction_reg == "rubin_ampuero") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); } else if (friction_reg == "simplified_prakash_clifton") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); } else { - AKANTU_DEBUG_ERROR("Do not know the following friction regularisation: " + AKANTU_ERROR("Do not know the following friction regularisation: " << friction_reg); } } // Friction Law: Linear Cohesive else if (friction_law == "linear_cohesive") { if (friction_reg == "no_regularisation") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); } else if (friction_reg == "rubin_ampuero") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); } else if (friction_reg == "simplified_prakash_clifton") { if (is_ntn_contact) friction = new NTNFriction(contact); else friction = new NTRFFriction(contact); } else { - AKANTU_DEBUG_ERROR("Do not know the following friction regularisation: " + AKANTU_ERROR("Do not know the following friction regularisation: " << friction_reg); } } else { - AKANTU_DEBUG_ERROR("Do not know the following friction law: " + AKANTU_ERROR("Do not know the following friction law: " << friction_law); } AKANTU_DEBUG_OUT(); return friction; } } // namespace akantu diff --git a/src/common/aka_array.cc b/src/common/aka_array.cc index a9950241f..d014ca3ab 100644 --- a/src/common/aka_array.cc +++ b/src/common/aka_array.cc @@ -1,124 +1,124 @@ /** * @file aka_array.cc * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Aug 18 2015 * * @brief Implementation of akantu::Array * * @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 #include /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_common.hh" namespace akantu { /* -------------------------------------------------------------------------- */ /* Functions ArrayBase */ /* -------------------------------------------------------------------------- */ ArrayBase::ArrayBase(ID id) : id(std::move(id)), allocated_size(0), size_(0), nb_component(1), size_of_type(0) {} /* -------------------------------------------------------------------------- */ ArrayBase::~ArrayBase() = default; /* -------------------------------------------------------------------------- */ void ArrayBase::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "ArrayBase [" << std::endl; stream << space << " + size : " << size_ << std::endl; stream << space << " + nb component : " << nb_component << std::endl; stream << space << " + allocated size : " << allocated_size << std::endl; Real mem_size = (allocated_size * nb_component * size_of_type) / 1024.; stream << space << " + size of type : " << size_of_type << "B" << std::endl; stream << space << " + memory allocated : " << mem_size << "kB" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ template <> UInt Array::find(const Real & elem) const { AKANTU_DEBUG_IN(); Real epsilon = std::numeric_limits::epsilon(); auto it = std::find_if(begin(), end(), [&elem, &epsilon](auto && a) { return std::abs(a - elem) <= epsilon; }); AKANTU_DEBUG_OUT(); return (it != end()) ? end() - it : UInt(-1); } /* -------------------------------------------------------------------------- */ template <> Array & Array::operator*=(__attribute__((unused)) const ElementType & alpha) { - AKANTU_DEBUG_TO_IMPLEMENT(); + AKANTU_TO_IMPLEMENT(); return *this; } template <> Array & Array:: operator-=(__attribute__((unused)) const Array & vect) { - AKANTU_DEBUG_TO_IMPLEMENT(); + AKANTU_TO_IMPLEMENT(); return *this; } template <> Array & Array:: operator+=(__attribute__((unused)) const Array & vect) { - AKANTU_DEBUG_TO_IMPLEMENT(); + AKANTU_TO_IMPLEMENT(); return *this; } template <> Array & Array::operator*=(__attribute__((unused)) const char & alpha) { - AKANTU_DEBUG_TO_IMPLEMENT(); + AKANTU_TO_IMPLEMENT(); return *this; } template <> Array & Array::operator-=(__attribute__((unused)) const Array & vect) { - AKANTU_DEBUG_TO_IMPLEMENT(); + AKANTU_TO_IMPLEMENT(); return *this; } template <> Array & Array::operator+=(__attribute__((unused)) const Array & vect) { - AKANTU_DEBUG_TO_IMPLEMENT(); + AKANTU_TO_IMPLEMENT(); return *this; } } // namespace akantu diff --git a/src/common/aka_array_tmpl.hh b/src/common/aka_array_tmpl.hh index 87082c357..5f779890e 100644 --- a/src/common/aka_array_tmpl.hh +++ b/src/common/aka_array_tmpl.hh @@ -1,1280 +1,1280 @@ /** * @file aka_array_tmpl.hh * * @author Nicolas Richart * * @date creation: Thu Jul 15 2010 * @date last modification: Fri Jan 22 2016 * * @brief Inline functions of the classes Array and ArrayBase * * @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 . * */ /* -------------------------------------------------------------------------- */ /* Inline Functions Array */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AKA_ARRAY_TMPL_HH__ #define __AKANTU_AKA_ARRAY_TMPL_HH__ namespace akantu { namespace debug { struct ArrayException : public Exception {}; } // namespace debug /* -------------------------------------------------------------------------- */ template inline T & Array::operator()(UInt i, UInt j) { AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty"); AKANTU_DEBUG_ASSERT((i < size_) && (j < nb_component), "The value at position [" << i << "," << j << "] is out of range in array \"" << id << "\""); return values[i * nb_component + j]; } /* -------------------------------------------------------------------------- */ template inline const T & Array::operator()(UInt i, UInt j) const { AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty"); AKANTU_DEBUG_ASSERT((i < size_) && (j < nb_component), "The value at position [" << i << "," << j << "] is out of range in array \"" << id << "\""); return values[i * nb_component + j]; } template inline T & Array::operator[](UInt i) { AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty"); AKANTU_DEBUG_ASSERT((i < size_ * nb_component), "The value at position [" << i << "] is out of range in array \"" << id << "\""); return values[i]; } /* -------------------------------------------------------------------------- */ template inline const T & Array::operator[](UInt i) const { AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty"); AKANTU_DEBUG_ASSERT((i < size_ * nb_component), "The value at position [" << i << "] is out of range in array \"" << id << "\""); return values[i]; } /* -------------------------------------------------------------------------- */ /** * append a tuple to the array with the value value for all components * @param value the new last tuple or the array will contain nb_component copies * of value */ template inline void Array::push_back(const T & value) { resizeUnitialized(size_ + 1, true, value); } /* -------------------------------------------------------------------------- */ /** * append a tuple to the array * @param new_elem a C-array containing the values to be copied to the end of * the array */ // template // inline void Array::push_back(const T new_elem[]) { // UInt pos = size_; // resizeUnitialized(size_ + 1, false); // T * tmp = values + nb_component * pos; // std::uninitialized_copy(new_elem, new_elem + nb_component, tmp); // } /* -------------------------------------------------------------------------- */ #ifndef SWIG /** * append a matrix or a vector to the array * @param new_elem a reference to a Matrix or Vector */ template template