diff --git a/examples/phase_field/phase_field_notch.cc b/examples/phase_field/phase_field_notch.cc index ff02b59d3..e6f8e0c65 100644 --- a/examples/phase_field/phase_field_notch.cc +++ b/examples/phase_field/phase_field_notch.cc @@ -1,93 +1,93 @@ /** * @file phase_field_static_2d.cc * * @author Mohit Pundir * * @date creation: Mon Oct 1 2018 * * @brief test of the class PhaseFieldModel on the 2d square * * @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 "non_linear_solver.hh" #include "phase_field_model.hh" #include "solid_mechanics_model.hh" #include "solid_phase_coupler.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; const UInt spatial_dimension = 2; /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { initialize("material_notch.dat", argc, argv); // create mesh Mesh mesh(spatial_dimension); mesh.read("square_notch.msh"); //Phase field model initialization PhaseFieldModel pfm(mesh); pfm.initFull(_analysis_method = _static); // solid mechanics model initialization SolidMechanicsModel smm(mesh); smm.initFull(_analysis_method = _static); smm.applyBC(BC::Dirichlet::FixedValue(0., _y), "bottom"); smm.applyBC(BC::Dirichlet::FixedValue(0., _x), "left"); - smm.setBaseName( "square_notch"); + smm.setBaseName( "square_notch_modified"); smm.addDumpFieldVector( "displacement"); smm.addDumpFieldVector( "internal_force"); smm.addDumpField( "stress"); smm.addDumpField( "grad_u"); smm.addDumpField( "damage"); smm.addDumpField( "blocked_dofs"); smm.dump(); auto & smm_solver = smm.getNonLinearSolver(); smm_solver.set("max_iterations", 1000); smm_solver.set("threshold", 1e-8); smm_solver.set("convergence_type", _scc_residual); // coupling of models SolidPhaseCoupler coupler(smm, pfm); - UInt nbSteps = 1500; + UInt nbSteps = 200; Real increment = 1.e-5; for (UInt s = 1; s < nbSteps; ++s) { smm.applyBC(BC::Dirichlet::IncrementValue(increment, _y), "top"); //smm.solveStep(); coupler.solve(); smm.dump(); std::cout << "Step " << s << "/" << nbSteps << std::endl; } finalize(); return EXIT_SUCCESS; } diff --git a/src/model/model_couplers/solid_phase_coupler.cc b/src/model/model_couplers/solid_phase_coupler.cc index bcd2b8597..018fa281c 100644 --- a/src/model/model_couplers/solid_phase_coupler.cc +++ b/src/model/model_couplers/solid_phase_coupler.cc @@ -1,209 +1,214 @@ /** * @file solid_phase_coupler.cc * * @author Mohit Pundir * * @date creation: Fri Sep 28 2018 * @date last modification: Fri Sep 28 2018 * * @brief class for coupling of solid mechancis and phasefield model * * @section LICENSE * * Copyright (©) 2010-2018 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_phase_coupler.hh" #include "non_linear_solver.hh" /* -------------------------------------------------------------------------- */ namespace akantu { template SolidPhaseCoupler::SolidPhaseCoupler(SolidType & solid, PhaseType & phase) : solid(solid), phase(phase) { this->spatial_dimension = solid.getMesh().getSpatialDimension(); } /* -------------------------------------------------------------------------- */ template SolidPhaseCoupler::~SolidPhaseCoupler() { } /* -------------------------------------------------------------------------- */ template void SolidPhaseCoupler::computeDamageOnQuadPoints(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); auto & fem = phase.getFEEngine(); auto & mesh = phase.getMesh(); switch (spatial_dimension) { case 1: { auto & mat = static_cast &>(solid.getMaterial(0)); auto & damage = mat.getDamage(); for (auto & type: mesh.elementTypes(this->spatial_dimension, ghost_type)) { auto & damage_on_qpoints_vect = damage(type, ghost_type); fem.interpolateOnIntegrationPoints(phase.getDamage(), damage_on_qpoints_vect, 1, type, ghost_type); } break; } case 2: { auto & mat = static_cast &>(solid.getMaterial(0)); auto & damage = mat.getDamage(); for (auto & type: mesh.elementTypes(this->spatial_dimension, ghost_type)) { auto & damage_on_qpoints_vect = damage(type, ghost_type); fem.interpolateOnIntegrationPoints(phase.getDamage(), damage_on_qpoints_vect, 1, type, ghost_type); } break; } default: auto & mat = static_cast &>(solid.getMaterial(0)); break; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void SolidPhaseCoupler::computeStrainOnQuadPoints(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); auto & mesh = solid.getMesh(); auto & strain_on_qpoints = phase.getStrain(); auto & gradu_on_qpoints = solid.getMaterial(0).getGradU(); for (auto & type: mesh.elementTypes(spatial_dimension, ghost_type)) { auto & strain_on_qpoints_vect = strain_on_qpoints(type, ghost_type); auto & gradu_on_qpoints_vect = gradu_on_qpoints(type, ghost_type); for (auto && values: zip(make_view(strain_on_qpoints_vect, this->spatial_dimension, this->spatial_dimension), make_view(gradu_on_qpoints_vect, this->spatial_dimension, this->spatial_dimension))) { auto & strain = std::get<0>(values); auto & grad_u = std::get<1>(values); this->gradUToEpsilon(grad_u, strain); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void SolidPhaseCoupler::solve() { this->convergence = true; UInt iter = 0; UInt max_iter = 10; while(iter < max_iter) { auto u_old = solid.getDisplacement(); auto d_old = phase.getDamage(); std::cerr << "---- solving solid model ------ \n"; solid.solveStep(); this->computeStrainOnQuadPoints(_not_ghost); std::cerr << "---- solving phasefield model ------ \n"; phase.solveStep(); this->computeDamageOnQuadPoints(_not_ghost); auto u_new = solid.getDisplacement(); auto d_new = phase.getDamage(); + //auto u_new = solid.getDOFManager().getDOFs("displacement"); + //auto u_old = solid.getDOFManager().getPreviousDOFs("displacement"); + + //auto d_new = phase.getDOFManager().getDOFs("damage"); + //auto d_old = phase.getDOFManager().getPreviousDOFs("damage"); - //this->testConvergence(u_new, u_old, d_new, d_old); + this->checkConvergence(u_new, u_old, d_new, d_old); if (this->convergence) { break; } iter++; } } /* -------------------------------------------------------------------------- */ template void SolidPhaseCoupler::gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon) { for (UInt i=0; i < this->spatial_dimension; ++i) { for (UInt j = 0; j < this->spatial_dimension; ++j) epsilon(i, j) = 0.5 * (grad_u(i, j) + grad_u(j, i)); } } /* -------------------------------------------------------------------------- */ template -bool SolidPhaseCoupler::testConvergence(Array & u_new, Array & u_old, Array & d_new, Array & d_old) { +bool SolidPhaseCoupler::checkConvergence(Array & u_new, Array & u_old, Array & d_new, Array & d_old) { const Array & blocked_dofs = solid.getBlockedDOFs(); UInt nb_degree_of_freedom = u_new.size(); auto u_n_it = u_new.begin(); auto u_o_it = u_old.begin(); auto bld_it = blocked_dofs.begin(); Real norm = 0; for (UInt n = 0; n < nb_degree_of_freedom; ++n, ++u_n_it, ++u_o_it, ++bld_it) { if ((!*bld_it)) { norm += (*u_n_it - *u_o_it) * (*u_n_it - *u_o_it); } } norm = std::sqrt(norm); auto d_n_it = d_new.begin(); auto d_o_it = d_old.begin(); nb_degree_of_freedom = d_new.size(); Real norm2 = 0; for (UInt i = 0; i < nb_degree_of_freedom; ++i) { norm2 += (*d_n_it - *d_o_it); } norm2 =std::sqrt(norm2); Real error = std::max(norm, norm2); Real tolerance = 1e-8; if (error < tolerance) { this->convergence = true; } } template class SolidPhaseCoupler; } // akantu diff --git a/src/model/model_couplers/solid_phase_coupler.hh b/src/model/model_couplers/solid_phase_coupler.hh index 193a9d56b..599e9c2f6 100644 --- a/src/model/model_couplers/solid_phase_coupler.hh +++ b/src/model/model_couplers/solid_phase_coupler.hh @@ -1,105 +1,105 @@ /** * @file solid_phase_coupler.hh * * @author Mohit Pundir * * @date creation: Fri Sep 28 2018 * @date last modification: Fri Sep 28 2018 * * @brief class for coupling of solid mechancis and phasefield model * * @section LICENSE * * Copyright (©) 2010-2018 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 "solid_mechanics_model.hh" #include "phase_field_model.hh" #include "material.hh" #include "material_phasefield.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_PHASE_COUPLER_HH__ #define __AKANTU_SOLID_PHASE_COUPLER_HH__ namespace akantu { template class SolidPhaseCoupler { /* ------------------------------------------------------------------------ */ /* Constructor/Destructors */ /* ------------------------------------------------------------------------ */ private: public: SolidPhaseCoupler(SolidType &, PhaseType &); ~SolidPhaseCoupler(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// computes damage on quad points for solid mechanics model from /// damage array from phasefield model void computeDamageOnQuadPoints(const GhostType &); /// computes strain on quadrature points for phasefield model from /// displacement gradient from solid mechanics model void computeStrainOnQuadPoints(const GhostType & ghost_type); /// solve the coupled model void solve(); private: /// computes small strain from displacement gradient void gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon); /// test the convergence criteria - bool testConvergence(Array &, Array &, Array &, Array &); + bool checkConvergence(Array &, Array &, Array &, Array &); /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ private: /// SolidType & solid; /// PhaseType & phase; /// Spatial dimension of models UInt spatial_dimension; /// bool convergence; }; } // akantu #endif /* __AKANTU_SOLID_PHASE_COUPLER_HH__ */ diff --git a/src/model/phase_field/phase_field_model.cc b/src/model/phase_field/phase_field_model.cc index d174216c4..9b8bb8dbb 100644 --- a/src/model/phase_field/phase_field_model.cc +++ b/src/model/phase_field/phase_field_model.cc @@ -1,893 +1,906 @@ /** * @file phase_field_model.cc * * @author Mohit Pundir * * @date creation: Wed Aug 01 2018 * @date last modification: Wed Aug 01 2018 * * @brief Implementation of PhaseFieldModel class * * @section LICENSE * * Copyright (©) 2010-2018 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 "phase_field_model.hh" #include "dumpable_inline_impl.hh" #include "element_synchronizer.hh" #include "generalized_trapezoidal.hh" #include "fe_engine_template.hh" #include "generalized_trapezoidal.hh" #include "group_manager_inline_impl.cc" #include "integrator_gauss.hh" #include "mesh.hh" #include "parser.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_element_partition.hh" #include "dumper_elemental_field.hh" #include "dumper_internal_material_field.hh" #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { namespace phasefield { namespace details { class ComputeDamageFunctor { public: ComputeDamageFunctor(const PhaseFieldModel & model) : model(model) {}; // material functor can be created void operator()(Matrix & rho, const Element &) const { rho.set(1.); } private: const PhaseFieldModel & model; }; } } /* -------------------------------------------------------------------------- */ PhaseFieldModel::PhaseFieldModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, ModelType::_phase_field_model, dim, id, memory_id), + BoundaryCondition(), damage_on_qpoints( "damage_on_qpoints", id), damage_energy_on_qpoints( "damage_energy_on_qpoints", id), damage_energy_density_on_qpoints( "damage_energy_density_on_qpoints", id), damage_gradient( "damage_gradient", id), strain_on_qpoints( "strain_on_qpoints", id), driving_force_on_qpoints( "driving_force_on_qpoints", id), phi_history_on_qpoints( "phi_history_on_qpoints", id){ AKANTU_DEBUG_IN(); this->initDOFManager(); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, _gst_pfm_damage); this->registerSynchronizer(synchronizer, _gst_pfm_gradient_damage); } this->registerFEEngineObject(id + ":fem", mesh, Model::spatial_dimension); #ifdef AKANTU_USE_IOHELPER this->mesh.registerDumper("phase_field", id, true); this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); #endif // AKANTU_USE_IOHELPER this->registerParam("l0", l_0 , 0., _pat_parsmod, "length scale"); this->registerParam("gc", g_c , _pat_parsmod, "critical local fracture energy density"); this->registerParam("E" , E , _pat_parsmod, "Young's modulus"); this->registerParam("nu", nu , _pat_parsmod, "Poisson ratio"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ PhaseFieldModel::~PhaseFieldModel() = default; /* -------------------------------------------------------------------------- */ MatrixType PhaseFieldModel::getMatrixType(const ID & matrix_id) { if (matrix_id == "K" or matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::initModel() { auto & fem = this->getFEEngine(); fem.initShapeFunctions(_not_ghost); fem.initShapeFunctions(_ghost); damage_on_qpoints.initialize(fem, _nb_component = 1); damage_energy_on_qpoints.initialize(fem, _nb_component = spatial_dimension * spatial_dimension); damage_energy_density_on_qpoints.initialize(fem, _nb_component = 1); damage_gradient.initialize(fem, _nb_component = spatial_dimension); strain_on_qpoints.initialize(fem, _nb_component = spatial_dimension * spatial_dimension); driving_force_on_qpoints.initialize(fem, _nb_component = 1); phi_history_on_qpoints.initialize(fem, _nb_component = 1); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); readMaterials(); + this->initBC(*this, *damage, *external_force); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::readMaterials() { auto sect = this->getParserSection(); if (not std::get<1>(sect)) { const auto & section = std::get<0>(sect); this->parseSection(section); } Matrix d(spatial_dimension, spatial_dimension); d.eye(g_c * l_0); damage_energy_on_qpoints.set(d); this->updateInternalParameters(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::updateInternalParameters() { this->lambda = this->nu * this->E / ((1 + this->nu) * (1 - 2 * this->nu)); this->mu = this->E / (2 * (1 + this->nu)); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleDamageMatrix(); this->assembleDamageGradMatrix(); } else { AKANTU_ERROR("Unknown Matrix ID for PhaseFieldModel : " << matrix_id); } } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::predictor() { //AKANTU_TO_IMPLEMENT(); } + +/* -------------------------------------------------------------------------- */ +void PhaseFieldModel::corrector() { + //AKANTU_TO_IMPLEMENT(); +} + + /* -------------------------------------------------------------------------- */ void PhaseFieldModel::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType) { DOFManager & dof_manager = this->getDOFManager(); this->allocNodalField(this->damage, 1, "damage"); this->allocNodalField(this->external_force, 1, "external_force"); this->allocNodalField(this->internal_force, 1, "internal_force"); this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs"); + this->allocNodalField(this->previous_damage, 1, "previous_damage"); + this->allocNodalField(this->damage_increment, 1, "damage_increment"); if (!dof_manager.hasDOFs("damage")) { dof_manager.registerDOFs("damage", *this->damage, _dst_nodal); dof_manager.registerBlockedDOFs("damage", *this->blocked_dofs); + dof_manager.registerDOFsIncrement("damage", + *this->damage_increment); + dof_manager.registerDOFsPrevious("damage", + *this->previous_damage); } if (time_step_solver_type == _tsst_dynamic) { AKANTU_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ FEEngine & PhaseFieldModel::getFEEngineBoundary(const ID & name) { return dynamic_cast(getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ std::tuple PhaseFieldModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", _tsst_dynamic_lumped); } case _static: { return std::make_tuple("static", _tsst_static); } case _implicit_dynamic: { return std::make_tuple("implicit", _tsst_dynamic); } default: return std::make_tuple("unknown", _tsst_not_defined); } } /* -------------------------------------------------------------------------- */ ModelSolverOptions PhaseFieldModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case _tsst_static: { options.non_linear_solver_type = _nls_linear; options.integration_scheme_type["damage"] = _ist_pseudo_time; options.solution_type["damage"] = IntegrationScheme::_not_defined; break; } case _tsst_dynamic: { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["damage"] = _ist_backward_euler; options.solution_type["damage"] = IntegrationScheme::_damage; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::beforeSolveStep() { this->computePhiHistoryOnQuadPoints(_not_ghost); this->computeDamageEnergyDensityOnQuadPoints(_not_ghost); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::afterSolveStep() { - - //for (auto & dam : *damage) { - // dam = std::min(1., 2*dam -dam * dam); - //} - + for (auto & dam : *damage) { + dam = std::min(1., 2*dam -dam * dam); + } } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleDamageMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new damage matrix"); if (!this->getDOFManager().hasMatrix("K")) { this->getDOFManager().getNewMatrix("K", getMatrixType("K")); } this->getDOFManager().clearMatrix("K"); switch (mesh.getSpatialDimension()) { case 1: this->assembleDamageMatrix<1>(_not_ghost); break; case 2: this->assembleDamageMatrix<2>(_not_ghost); break; case 3: this->assembleDamageMatrix<3>(_not_ghost); break; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleDamageGradMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new damage gradient matrix"); if (!this->getDOFManager().hasMatrix("K")) { this->getDOFManager().getNewMatrix("K", getMatrixType("K")); } switch (mesh.getSpatialDimension()) { case 1: this->assembleDamageGradMatrix<1>(_not_ghost); break; case 2: this->assembleDamageGradMatrix<2>(_not_ghost); break; case 3: this->assembleDamageGradMatrix<3>(_not_ghost); break; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void PhaseFieldModel::assembleDamageMatrix(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); auto & fem = getFEEngineClass(); for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type)) { auto nb_element = mesh.getNbElement(type, ghost_type); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); auto nt_b_n = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "N^t*b*N"); auto & damage_energy_density_on_qpoints_vect = damage_energy_density_on_qpoints(type, ghost_type); // damage_energy_on_qpoints = gc/l0 + phi = scalar fem.computeNtbN(damage_energy_density_on_qpoints_vect, *nt_b_n, 2, type, ghost_type); /// compute @f$ K_{\grad d} = \int_e \mathbf{N}^t * \mathbf{w} * /// \mathbf{N}@f$ auto K_n = std::make_unique>( nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_n"); fem.integrate(*nt_b_n, *K_n, nb_nodes_per_element * nb_nodes_per_element, type, ghost_type); this->getDOFManager().assembleElementalMatricesToMatrix( "K", "damage", *K_n, type, ghost_type, _symmetric); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void PhaseFieldModel::assembleDamageGradMatrix(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); auto & fem = this->getFEEngine(); for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type)) { auto nb_element = mesh.getNbElement(type, ghost_type); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); auto bt_d_b = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "B^t*D*B"); auto & damage_energy_on_qpoints_vect = damage_energy_on_qpoints(type, ghost_type); // damage_energy_on_qpoints = gc*l0 = scalar fem.computeBtDB(damage_energy_on_qpoints_vect, *bt_d_b, 2, type, ghost_type); /// compute @f$ K_{\grad d} = \int_e \mathbf{B}^t * \mathbf{W} * /// \mathbf{B}@f$ auto K_b = std::make_unique>( nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_b"); fem.integrate(*bt_d_b, *K_b, nb_nodes_per_element * nb_nodes_per_element, type, ghost_type); this->getDOFManager().assembleElementalMatricesToMatrix( "K", "damage", *K_b, type, ghost_type, _symmetric); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::computeDamageEnergyDensityOnQuadPoints( const GhostType & ghost_type) { AKANTU_DEBUG_IN(); for (auto & type: mesh.elementTypes(spatial_dimension, ghost_type)) { for (auto && values : zip(make_view(damage_energy_density_on_qpoints(type, ghost_type)), make_view(phi_history_on_qpoints(type, ghost_type)))) { auto & dam_energy_density = std::get<0>(values); auto & phi_history = std::get<1>(values); dam_energy_density = g_c/l_0 + 2.0 * phi_history; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::computePhiHistoryOnQuadPoints( const GhostType & ghost_type) { Matrix strain_plus( spatial_dimension, spatial_dimension); Matrix strain_minus( spatial_dimension, spatial_dimension); Matrix strain_dir( spatial_dimension, spatial_dimension); Matrix strain_diag_plus( spatial_dimension, spatial_dimension); Matrix strain_diag_minus(spatial_dimension, spatial_dimension); Vector strain_values(spatial_dimension); Real trace_plus, trace_minus, phi_plus; for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) { for (auto && values : zip(make_view(strain_on_qpoints(type, ghost_type), spatial_dimension, spatial_dimension), make_view(phi_history_on_qpoints(type, ghost_type)))) { auto & strain = std::get<0>(values); auto & phi_history = std::get<1>(values); strain_plus.clear(); strain_minus.clear(); strain_dir.clear(); strain_values.clear(); strain_diag_plus.clear(); strain_diag_minus.clear(); strain.eig(strain_values, strain_dir); for (UInt i=0; i < spatial_dimension; i++) { strain_diag_plus(i, i) = std::max(Real(0.), strain_values(i)); strain_diag_minus(i, i) = std::min(Real(0.), strain_values(i)); } Matrix mat_tmp( spatial_dimension, spatial_dimension); Matrix sigma_plus( spatial_dimension, spatial_dimension); Matrix sigma_minus(spatial_dimension, spatial_dimension); - mat_tmp.mul(strain_diag_plus, strain_dir); + mat_tmp.mul(strain_diag_plus, strain_dir); strain_plus.mul(strain_dir, mat_tmp); mat_tmp.mul(strain_diag_minus, strain_dir); strain_minus.mul(strain_dir, mat_tmp); trace_plus = std::max(Real(0.), strain.trace()); trace_minus = std::min(Real(0.), strain.trace()); for (UInt i=0; i < spatial_dimension; i++) { for (UInt j=0; j < spatial_dimension; j++) { sigma_plus(i, j) = (i==j) * lambda * trace_plus + 2 * mu * strain_plus(i, j); sigma_minus(i, j) = (i==j) * lambda * trace_minus + 2 * mu * strain_minus(i, j); } } phi_plus = 1./2 * sigma_plus.doubleDot(strain); if (phi_plus > phi_history) { phi_history = phi_plus; } } } } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleResidual() { AKANTU_DEBUG_IN(); this->assembleInternalForces(); this->getDOFManager().assembleToResidual("damage", *this->external_force, 1); this->getDOFManager().assembleToResidual("damage", *this->internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleInternalForces() { AKANTU_DEBUG_IN(); this->internal_force->clear(); this->synchronize(_gst_pfm_damage); auto & fem = this->getFEEngine(); for (auto ghost_type: ghost_types) { computeDrivingForce(ghost_type); for (auto type: mesh.elementTypes(spatial_dimension, ghost_type)) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto & driving_force_on_qpoints_vect = driving_force_on_qpoints(type, ghost_type); UInt nb_quad_points = driving_force_on_qpoints_vect.size(); Array nt_driving_force(nb_quad_points, nb_nodes_per_element); fem.computeNtb(driving_force_on_qpoints_vect, nt_driving_force, type, ghost_type); UInt nb_elements = mesh.getNbElement(type, ghost_type); Array int_nt_driving_force(nb_elements, nb_nodes_per_element); fem.integrate(nt_driving_force, int_nt_driving_force, nb_nodes_per_element, type, ghost_type); this->getDOFManager().assembleElementalArrayLocalArray( int_nt_driving_force, *this->internal_force, type, ghost_type, 1); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleLumpedMatrix(const ID & matrix_id) { } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper("phase_field").setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::computeDrivingForce(const GhostType & ghost_type) { for (auto & type: mesh.elementTypes(spatial_dimension, ghost_type)) { for (auto && values: zip(make_view(phi_history_on_qpoints(type, ghost_type)), make_view(driving_force_on_qpoints(type, ghost_type)) )) { auto & phi_history = std::get<0>(values); auto & driving_force = std::get<1>(values); driving_force = 2.0 * phi_history; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::computeDamageOnQuadPoints(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); auto & fem = this->getFEEngine(); for (auto & type: mesh.elementTypes(spatial_dimension, ghost_type)) { auto & damage_on_qpoints_vect = damage_on_qpoints(type, ghost_type); fem.interpolateOnIntegrationPoints(*damage, damage_on_qpoints_vect, 1, type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt PhaseFieldModel::getNbData(const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = 0; for (const Element & el : elements) { nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } switch (tag) { case _gst_material_id: { size += elements.size() * sizeof(UInt); break; } case _gst_pfm_damage: { size += nb_nodes_per_element * sizeof(Real); // damage break; } case _gst_pfm_gradient_damage: { size += getNbIntegrationPoints(elements) * spatial_dimension * sizeof(Real); size += nb_nodes_per_element * sizeof(Real); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { switch (tag) { case _gst_pfm_damage: { packNodalDataHelper(*damage, buffer, elements, mesh); break; } case _gst_pfm_gradient_damage: { packElementalDataHelper(damage_gradient, buffer, elements, true, getFEEngine()); packNodalDataHelper(*damage, buffer, elements, mesh); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { switch (tag) { case _gst_pfm_damage: { unpackNodalDataHelper(*damage, buffer, elements, mesh); break; } case _gst_pfm_gradient_damage: { unpackElementalDataHelper(damage_gradient, buffer, elements, true, getFEEngine()); unpackNodalDataHelper(*damage, buffer, elements, mesh); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } /* -------------------------------------------------------------------------- */ UInt PhaseFieldModel::getNbData(const Array & indexes, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes = indexes.size(); switch (tag) { case _gst_pfm_damage: { size += nb_nodes * sizeof(Real); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::packData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); for (auto index : indexes) { switch (tag) { case _gst_pfm_damage: { buffer << (*damage)(index); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::unpackData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); for (auto index : indexes) { switch (tag) { case _gst_pfm_damage: { buffer >> (*damage)(index); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER dumper::Field * PhaseFieldModel::createNodalFieldBool( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag ) { std::map *> uint_nodal_fields; uint_nodal_fields["blocked_dofs"] = blocked_dofs; dumper::Field * field = nullptr; field = mesh.createNodalField(uint_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ dumper::Field * PhaseFieldModel::createNodalFieldReal( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { if (field_name == "capacity_lumped") { AKANTU_EXCEPTION("Capacity lumped is a nodal field now stored in the DOF manager." "Therefore it cannot be used by a dumper anymore"); } std::map *> real_nodal_fields; real_nodal_fields["damage"] = damage; real_nodal_fields["external_force"] = external_force; real_nodal_fields["internal_force"] = internal_force; dumper::Field * field = mesh.createNodalField(real_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ dumper::Field * PhaseFieldModel::createElementalField( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag, __attribute__((unused)) const UInt & spatial_dimension, const ElementKind & element_kind) { dumper::Field * field = nullptr; if (field_name == "partitions") field = mesh.createElementalField( mesh.getConnectivities(), group_name, this->spatial_dimension, element_kind); else if (field_name == "damage_gradient") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(damage_gradient, element_kind); field = mesh.createElementalField( damage_gradient, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } else if (field_name == "damage_energy") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(damage_energy_on_qpoints, element_kind); field = mesh.createElementalField( damage_energy_on_qpoints, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } return field; } /* -------------------------------------------------------------------------- */ #else /* -------------------------------------------------------------------------- */ dumper::Field * PhaseFieldModel::createElementalField( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag, __attribute__((unused)) const ElementKind & element_kind) { return nullptr; } /* -------------------------------------------------------------------------- */ dumper::Field * PhaseFieldModel::createNodalFieldBool( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return nullptr; } /* -------------------------------------------------------------------------- */ dumper::Field * PhaseFieldModel::createNodalFieldReal( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void PhaseFieldModel::dump(const std::string & dumper_name) { mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::dump(const std::string & dumper_name, UInt step) { mesh.dump(dumper_name, step); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::dump(const std::string & dumper_name, Real time, UInt step) { mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::dump() { mesh.dump(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::dump(UInt step) { mesh.dump(step); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::dump(Real time, UInt step) { mesh.dump(time, step); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "Phase Field Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << Model::spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; damage->printself(stream, indent + 2); external_force->printself(stream, indent + 2); internal_force->printself(stream, indent + 2); blocked_dofs->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + material information [" << std::endl; stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + materials [" << std::endl; stream << space << " length scale parameter : " << l_0 << std::endl; stream << space << " critical energy release rate : " << g_c << std::endl; stream << space << " Young's modulus : " << E << std::endl; stream << space << " Poisson's ratio : " << nu << std::endl; stream << space << " Lame's first parameter : " << lambda << std::endl; stream << space << " Lame's second parameter : " << mu << std::endl; stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } } // akantu diff --git a/src/model/phase_field/phase_field_model.hh b/src/model/phase_field/phase_field_model.hh index ec652663c..8719a30fd 100644 --- a/src/model/phase_field/phase_field_model.hh +++ b/src/model/phase_field/phase_field_model.hh @@ -1,294 +1,315 @@ /** * @file phase_field_model.hh * * @author Mohit Pundir * * @date creation: Sun Jul 30 2018 * @date last modification: Mon Feb 05 2018 * * @brief Model of Phase Field * * @section LICENSE * * Copyright (©) 2010-2018 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 "boundary_condition.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PHASE_FIELD_MODEL_HH__ #define __AKANTU_PHASE_FIELD_MODEL_HH__ namespace akantu { template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu namespace akantu { class PhaseFieldModel : public Model, public DataAccessor, - public DataAccessor { + public DataAccessor, + public BoundaryCondition { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using FEEngineType = FEEngineTemplate; PhaseFieldModel(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "phase_field_model", const MemoryID & memory_id = 0); virtual ~PhaseFieldModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// generic function to initialize everything ready for explicit dynamics void initFullImpl(const ModelOptions & options) override; /// read one material file to instantiate all the materials void readMaterials(); /// allocate all vectors void initSolver(TimeStepSolverType, NonLinearSolverType) override; /// initialize the model void initModel() override; + /// predictor void predictor() override; + /// corrector + void corrector() override; + /// compute the heat flux void assembleResidual() override; /// get the type of matrix needed MatrixType getMatrixType(const ID &) override; /// callback to assemble a Matrix void assembleMatrix(const ID &) override; /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID &) override; std::tuple getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const; /// compute the internal forces void assembleInternalForces(); /// void updateInternalParameters(); /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /* ------------------------------------------------------------------------ */ /* Methods for static */ /* ------------------------------------------------------------------------ */ public: /// assemble damage matrix void assembleDamageMatrix(); /// assemble damage gradient matrix void assembleDamageGradMatrix(); /// coupling parameters damage and strains from solid mechanics model void setCouplingParameters(ElementTypeMapArray & strain_on_qpoints, Array & damage); private: /// assemble the damage matrix (w/ ghost type) template void assembleDamageMatrix(const GhostType & ghost_type); /// assemble the damage grad matrix (w/ ghost type) template void assembleDamageGradMatrix(const GhostType & ghost_type); private: /// compute vector strain history field for each quadrature point void computePhiHistoryOnQuadPoints(const GhostType & ghost_type); /// compute vector strain history field for each quadrature point void computeDamageEnergyDensityOnQuadPoints(const GhostType & ghost_type); /// compute driving force for each quadrature point void computeDrivingForce(const GhostType & ghost_type); /// compute the damage on quadrature points void computeDamageOnQuadPoints(const GhostType & ghost_type); /// compute the fracture energy Real computeFractureEnergyByNode(); /* ------------------------------------------------------------------------ */ /* Methods for static */ /* ------------------------------------------------------------------------ */ public: /// set the stable timestep void setTimeStep(Real time_step, const ID & solver_id="") override; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep() override; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & indexes, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: + /// return the dimension of the system space + AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); + + /// return the damage array AKANTU_GET_MACRO(Damage, *damage, Array &); + + /// AKANTU_GET_MACRO_NOT_CONST(Strain, strain_on_qpoints, ElementTypeMapArray &); - + + /// get the PhaseFieldModel::blocked_dofs vector + AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); + /* ------------------------------------------------------------------------ */ /* Dumpable Interface */ /* ------------------------------------------------------------------------ */ public: dumper::Field * createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) override; virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); protected: /* ------------------------------------------------------------------------ */ FEEngine & getFEEngineBoundary(const ID & name = "") override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: ///number of iterations UInt n_iter; /// damage array Array * damage{nullptr}; + /// damage array at the previous time step (used in finite deformation) + Array * previous_damage{nullptr}; + + /// increment of damage + Array * damage_increment{nullptr}; + /// damage field on quadrature points ElementTypeMapArray damage_on_qpoints; /// critical local damage energy on quadrature points for \mathbf{B}^t * \mathbf{W} * /// \mathbf{B}@f$ ElementTypeMapArray damage_energy_on_qpoints; /// critical local damage energy density on quadrature points for /// \mathbf{N}^t * \mathbf{w} * \mathbf{N}@f$ ElementTypeMapArray damage_energy_density_on_qpoints; /// the speed of change in damage ElementTypeMapArray damage_gradient; /// strain on quadrature points ElementTypeMapArray strain_on_qpoints; /// driving force on quadrature points for internal forces ElementTypeMapArray driving_force_on_qpoints; /// vector \phi plus on quadrature points ElementTypeMapArray phi_history_on_qpoints; /// boundary vector Array * blocked_dofs{nullptr}; /// external force vector Array * external_force{nullptr}; /// residuals array Array * internal_force{nullptr}; /// lengthscale parameter Real l_0; /// critical energy release rate Real g_c; /// Young's modulus Real E; /// Poisson ratio Real nu; /// Lame's first parameter Real lambda; /// Lame's second paramter Real mu; }; } // akantu /* ------------------------------------------------------------------------ */ /* inline functions */ /* ------------------------------------------------------------------------ */ #include "phase_field_model_inline_impl.cc" #endif diff --git a/test/test_model/CMakeLists.txt b/test/test_model/CMakeLists.txt index 2d7e0230f..9f3b49ffb 100644 --- a/test/test_model/CMakeLists.txt +++ b/test/test_model/CMakeLists.txt @@ -1,30 +1,31 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux # # @date creation: Fri Sep 03 2010 # @date last modification: Tue Jan 30 2018 # # @brief configuration for model tests # # @section LICENSE # # Copyright (©) 2010-2018 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 . # # @section DESCRIPTION # #=============================================================================== add_subdirectory(patch_tests) add_akantu_test(test_model_solver "Test for the solvers") add_akantu_test(test_solid_mechanics_model "Test for the solid mechanics model") add_akantu_test(test_structural_mechanics_model "Test for the structural mechanics model") add_akantu_test(test_non_local_toolbox "Test of the functionalities in the non-local toolbox") add_akantu_test(test_heat_transfer_model "Test heat transfer model") +add_akantu_test(test_phase_field_model "Test phase field model") diff --git a/test/test_model/CMakeLists.txt b/test/test_model/test_phase_field_model/CMakeLists.txt similarity index 58% copy from test/test_model/CMakeLists.txt copy to test/test_model/test_phase_field_model/CMakeLists.txt index 2d7e0230f..f6d80939e 100644 --- a/test/test_model/CMakeLists.txt +++ b/test/test_model/test_phase_field_model/CMakeLists.txt @@ -1,30 +1,43 @@ #=============================================================================== # @file CMakeLists.txt # -# @author Guillaume Anciaux +# @author Mohit Pundir # -# @date creation: Fri Sep 03 2010 -# @date last modification: Tue Jan 30 2018 +# @date creation: Thu Dec 20 2018 +# @date last modification: Thu Dec 20 2018 # -# @brief configuration for model tests +# @brief test for the phase field model # # @section LICENSE # # Copyright (©) 2010-2018 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 . # -# @section DESCRIPTION -# #=============================================================================== -add_subdirectory(patch_tests) -add_akantu_test(test_model_solver "Test for the solvers") -add_akantu_test(test_solid_mechanics_model "Test for the solid mechanics model") -add_akantu_test(test_structural_mechanics_model "Test for the structural mechanics model") -add_akantu_test(test_non_local_toolbox "Test of the functionalities in the non-local toolbox") -add_akantu_test(test_heat_transfer_model "Test heat transfer model") +add_mesh(test_one_element test_one_element.geo 2 1) +register_test(test_phase_field_model_2d + SOURCES test_phase_field_model_2d.cc + DEPENDS test_one_element + FILES_TO_COPY material.dat + PACKAGE phase_field + ) + +register_test(test_solid_one_element + SOURCES test_solid_one_element.cc + DEPENDS test_one_element + FILES_TO_COPY material_solid.dat + PACKAGE core + ) + +register_test(test_phase_solid_coupling + SOURCES test_phase_solid_coupling.cc + DEPENDS test_one_element + FILES_TO_COPY material_coupling.dat + PACKAGE phase_field core + ) diff --git a/test/test_model/test_phase_field_model/material.dat b/test/test_model/test_phase_field_model/material.dat new file mode 100644 index 000000000..98537c081 --- /dev/null +++ b/test/test_model/test_phase_field_model/material.dat @@ -0,0 +1,22 @@ +material elastic [ + name = plate + rho = 1. + E = 210.0 + nu = 0.3 + Plane_Stress = false +] + +material phasefield [ + name = plate + rho = 1. + E = 210.0 + nu = 0.3 + Plane_Stress = false +] + +model phase_field_model [ + l0 = 0.1 + gc = 5.0e-3 + E = 210.0 + nu = 0.3 +] \ No newline at end of file diff --git a/test/test_model/test_phase_field_model/material_coupling.dat b/test/test_model/test_phase_field_model/material_coupling.dat new file mode 100644 index 000000000..9eb1a0915 --- /dev/null +++ b/test/test_model/test_phase_field_model/material_coupling.dat @@ -0,0 +1,14 @@ +material phasefield [ + name = plate + rho = 1. + E = 210.0 + nu = 0.3 + Plane_Stress = false +] + +model phase_field_model [ + l0 = 0.1 + gc = 5.0e-3 + E = 210.0 + nu = 0.3 +] \ No newline at end of file diff --git a/test/test_model/test_phase_field_model/material_solid.dat b/test/test_model/test_phase_field_model/material_solid.dat new file mode 100644 index 000000000..8f185ee14 --- /dev/null +++ b/test/test_model/test_phase_field_model/material_solid.dat @@ -0,0 +1,7 @@ +material elastic [ + name = plate + rho = 1. + E = 210.0 + nu = 0.3 + Plane_Stress = false +] \ No newline at end of file diff --git a/test/test_model/test_phase_field_model/test_one_element.geo b/test/test_model/test_phase_field_model/test_one_element.geo new file mode 100644 index 000000000..840b30ca3 --- /dev/null +++ b/test/test_model/test_phase_field_model/test_one_element.geo @@ -0,0 +1,17 @@ +// Gmsh project created on Fri Dec 28 23:39:32 2018 +Point(1) = {-0.5, -0.5, 0, 100.0}; +Point(2) = {0.5, -0.5, 0, 100.0}; +Point(3) = {0.5, 0.5, 0, 100.0}; +Point(4) = {-0.5, 0.5, 0, 100.0}; + +Line(1) = {1, 2}; +Line(2) = {2, 3}; +Line(3) = {3, 4}; +Line(4) = {4, 1}; + +Line Loop(1) = {3, 4, 1, 2}; + +Plane Surface(1) = {1}; +Physical Surface("plate") = {1}; + +Recombine Surface{1}; \ No newline at end of file diff --git a/examples/phase_field/phase_field_notch.cc b/test/test_model/test_phase_field_model/test_phase_field_model_2d.cc similarity index 50% copy from examples/phase_field/phase_field_notch.cc copy to test/test_model/test_phase_field_model/test_phase_field_model_2d.cc index ff02b59d3..00664b470 100644 --- a/examples/phase_field/phase_field_notch.cc +++ b/test/test_model/test_phase_field_model/test_phase_field_model_2d.cc @@ -1,93 +1,98 @@ /** - * @file phase_field_static_2d.cc + * @file tets_phase_field_2d.cc * * @author Mohit Pundir * * @date creation: Mon Oct 1 2018 * * @brief test of the class PhaseFieldModel on the 2d square * * @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 "non_linear_solver.hh" #include "phase_field_model.hh" #include "solid_mechanics_model.hh" -#include "solid_phase_coupler.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; const UInt spatial_dimension = 2; /* -------------------------------------------------------------------------- */ +void applyStrainOnQuadPoints(PhaseFieldModel &); +/* -------------------------------------------------------------------------- */ -int main(int argc, char *argv[]) -{ - initialize("material_notch.dat", argc, argv); +int main(int argc, char *argv[]) { - // create mesh + initialize("material.dat", argc, argv); + Mesh mesh(spatial_dimension); - mesh.read("square_notch.msh"); - - //Phase field model initialization - PhaseFieldModel pfm(mesh); - pfm.initFull(_analysis_method = _static); - - // solid mechanics model initialization - SolidMechanicsModel smm(mesh); - smm.initFull(_analysis_method = _static); - - smm.applyBC(BC::Dirichlet::FixedValue(0., _y), "bottom"); - smm.applyBC(BC::Dirichlet::FixedValue(0., _x), "left"); - - smm.setBaseName( "square_notch"); - smm.addDumpFieldVector( "displacement"); - smm.addDumpFieldVector( "internal_force"); - smm.addDumpField( "stress"); - smm.addDumpField( "grad_u"); - smm.addDumpField( "damage"); - smm.addDumpField( "blocked_dofs"); - smm.dump(); - - auto & smm_solver = smm.getNonLinearSolver(); - smm_solver.set("max_iterations", 1000); - smm_solver.set("threshold", 1e-8); - smm_solver.set("convergence_type", _scc_residual); + mesh.read("test_one_element.msh"); + + PhaseFieldModel model(mesh); + model.initFull(_analysis_method = _static); + + model.setBaseName("one_element"); + model.addDumpFieldVector("damage"); + model.dump(); + + UInt nbSteps = 1000; + + for (UInt s = 0; s < nbSteps; ++s) { + applyStrainOnQuadPoints(model); + model.solveStep(); + model.dump(); + } - // coupling of models - SolidPhaseCoupler coupler(smm, pfm); + finalize(); + + return EXIT_SUCCESS; + +} + + +/* -------------------------------------------------------------------------- */ +void applyStrainOnQuadPoints(PhaseFieldModel & model) { + auto & strain_on_qpoints = model.getStrain(); - UInt nbSteps = 1500; - Real increment = 1.e-5; + auto & mesh = model.getMesh(); + auto & fem = model.getFEEngine(); - for (UInt s = 1; s < nbSteps; ++s) { - smm.applyBC(BC::Dirichlet::IncrementValue(increment, _y), "top"); - //smm.solveStep(); - coupler.solve(); - smm.dump(); - std::cout << "Step " << s << "/" << nbSteps << std::endl; + for (auto & type: mesh.elementTypes(spatial_dimension, _not_ghost) ) { + auto & strain_on_qpoints_vect = strain_on_qpoints(type, _not_ghost); + UInt nb_quad_points = fem.getNbIntegrationPoints(type, _not_ghost); + UInt nb_elements = mesh.getNbElement(type); - } + Array quad_coords(nb_elements * nb_quad_points, spatial_dimension); + fem.interpolateOnIntegrationPoints(mesh.getNodes(), quad_coords, spatial_dimension, type, _not_ghost); - finalize(); - return EXIT_SUCCESS; + for (auto && values: + zip(make_view(quad_coords), + make_view(strain_on_qpoints_vect, spatial_dimension, spatial_dimension))) { + auto & coord = std::get<0>(values); + auto & strain = std::get<1>(values); + + strain(1,1) += 1.e-4; + } + } } diff --git a/test/test_model/test_phase_field_model/test_phase_solid_coupling.cc b/test/test_model/test_phase_field_model/test_phase_solid_coupling.cc new file mode 100644 index 000000000..09f94f597 --- /dev/null +++ b/test/test_model/test_phase_field_model/test_phase_solid_coupling.cc @@ -0,0 +1,196 @@ +/** + * @file tets_phase_field_2d.cc + * + * @author Mohit Pundir + * + * @date creation: Mon Oct 1 2018 + * + * @brief test of the class PhaseFieldModel on the 2d square + * + * @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 "aka_common.hh" +#include "non_linear_solver.hh" +#include "solid_mechanics_model.hh" +#include "phase_field_model.hh" +#include "material.hh" +#include "material_phasefield.hh" +/* -------------------------------------------------------------------------- */ +#include +#include +/* -------------------------------------------------------------------------- */ + +using namespace akantu; +const UInt spatial_dimension = 2; + +/* -------------------------------------------------------------------------- */ +void applyDisplacement(SolidMechanicsModel &); +void computeStrainOnQuadPoints(SolidMechanicsModel &, PhaseFieldModel &, const GhostType &); +void computeDamageOnQuadPoints(SolidMechanicsModel &, PhaseFieldModel &, const GhostType &); +void gradUToEpsilon(const Matrix &, Matrix &); +/* -------------------------------------------------------------------------- */ + +int main(int argc, char *argv[]) { + + std::ofstream os("data.csv"); + os << "#disp stress damage" << std::endl; + + initialize("material_coupling.dat", argc, argv); + + Mesh mesh(spatial_dimension); + mesh.read("test_one_element.msh"); + + SolidMechanicsModel model(mesh); + model.initFull(_analysis_method = _static); + + PhaseFieldModel phase(mesh); + phase.initFull(_analysis_method = _static); + + model.setBaseName("phase_solid"); + model.addDumpField("stress"); + model.addDumpField("grad_u"); + model.addDumpField("damage"); + model.addDumpField("displacement"); + model.dump(); + + UInt nbSteps = 1000; + Real increment = 1e-4; + + auto & damage = model.getMaterial(0).getArray("damage", _quadrangle_4); + auto & stress = model.getMaterial(0).getArray("stress", _quadrangle_4); + + for (UInt s = 1; s < nbSteps; ++s) { + applyDisplacement(model); + + model.solveStep(); + computeStrainOnQuadPoints(model, phase, _not_ghost); + + phase.solveStep(); + computeDamageOnQuadPoints(model, phase, _not_ghost); + + os << s*increment << " " << stress(0, 0) << " " << damage(0) << std::endl; + + model.dump(); + } + + os.close(); + finalize(); + + return EXIT_SUCCESS; + +} + + +/* -------------------------------------------------------------------------- */ +void applyDisplacement(SolidMechanicsModel & model) { + auto & displacement = model.getDisplacement(); + + auto & positions = model.getMesh().getNodes(); + auto & blocked_dofs = model.getBlockedDOFs(); + + + for (UInt n = 0; n < model.getMesh().getNbNodes(); ++n) { + if (positions(n, 1) == -0.5) { + displacement(n, 0) = 0; + displacement(n, 1) = 0; + blocked_dofs(n, 0) = true; + blocked_dofs(n ,1) = true; + } + else { + displacement(n, 0) = 0; + displacement(n, 1) += 1.e-4; + blocked_dofs(n, 0) = true; + blocked_dofs(n ,1) = true; + } + } +} + +/* -------------------------------------------------------------------------- */ +void computeStrainOnQuadPoints(SolidMechanicsModel & solid, PhaseFieldModel & phase, const GhostType & ghost_type) { + + auto & mesh = solid.getMesh(); + + auto & strain_on_qpoints = phase.getStrain(); + auto & gradu_on_qpoints = solid.getMaterial(0).getGradU(); + + for (auto & type: mesh.elementTypes(spatial_dimension, ghost_type)) { + auto & strain_on_qpoints_vect = strain_on_qpoints(type, ghost_type); + auto & gradu_on_qpoints_vect = gradu_on_qpoints(type, ghost_type); + for (auto && values: + zip(make_view(strain_on_qpoints_vect, spatial_dimension, spatial_dimension), + make_view(gradu_on_qpoints_vect, spatial_dimension, spatial_dimension))) { + auto & strain = std::get<0>(values); + auto & grad_u = std::get<1>(values); + gradUToEpsilon(grad_u, strain); + } + + + } + +} + + +/* -------------------------------------------------------------------------- */ +void computeDamageOnQuadPoints(SolidMechanicsModel & solid, PhaseFieldModel & phase, const GhostType & ghost_type) { + + auto & fem = phase.getFEEngine(); + auto & mesh = phase.getMesh(); + + switch (spatial_dimension) { + case 1: { + auto & mat = static_cast &>(solid.getMaterial(0)); + auto & damage = mat.getDamage(); + for (auto & type: mesh.elementTypes(spatial_dimension, ghost_type)) { + auto & damage_on_qpoints_vect = damage(type, ghost_type); + fem.interpolateOnIntegrationPoints(phase.getDamage(), damage_on_qpoints_vect, + 1, type, ghost_type); + } + + break; + } + case 2: { + auto & mat = static_cast &>(solid.getMaterial(0)); + auto & damage = mat.getDamage(); + + for (auto & type: mesh.elementTypes(spatial_dimension, ghost_type)) { + auto & damage_on_qpoints_vect = damage(type, ghost_type); + fem.interpolateOnIntegrationPoints(phase.getDamage(), damage_on_qpoints_vect, + 1, type, ghost_type); + } + break; + } + default: + auto & mat = static_cast &>(solid.getMaterial(0)); + break; + } + +} + + +/* -------------------------------------------------------------------------- */ +void gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon) { + for (UInt i=0; i < spatial_dimension; ++i) { + for (UInt j = 0; j < spatial_dimension; ++j) + epsilon(i, j) = 0.5 * (grad_u(i, j) + grad_u(j, i)); + } +} diff --git a/examples/phase_field/phase_field_notch.cc b/test/test_model/test_phase_field_model/test_solid_one_element.cc similarity index 51% copy from examples/phase_field/phase_field_notch.cc copy to test/test_model/test_phase_field_model/test_solid_one_element.cc index ff02b59d3..e7aee5f20 100644 --- a/examples/phase_field/phase_field_notch.cc +++ b/test/test_model/test_phase_field_model/test_solid_one_element.cc @@ -1,93 +1,97 @@ /** - * @file phase_field_static_2d.cc + * @file tets_phase_field_2d.cc * * @author Mohit Pundir * * @date creation: Mon Oct 1 2018 * * @brief test of the class PhaseFieldModel on the 2d square * * @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 "non_linear_solver.hh" -#include "phase_field_model.hh" #include "solid_mechanics_model.hh" -#include "solid_phase_coupler.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; const UInt spatial_dimension = 2; /* -------------------------------------------------------------------------- */ +void applyDisplacement(SolidMechanicsModel &); +/* -------------------------------------------------------------------------- */ -int main(int argc, char *argv[]) -{ - initialize("material_notch.dat", argc, argv); +int main(int argc, char *argv[]) { - // create mesh + initialize("material_solid.dat", argc, argv); + Mesh mesh(spatial_dimension); - mesh.read("square_notch.msh"); - - //Phase field model initialization - PhaseFieldModel pfm(mesh); - pfm.initFull(_analysis_method = _static); - - // solid mechanics model initialization - SolidMechanicsModel smm(mesh); - smm.initFull(_analysis_method = _static); + mesh.read("test_one_element.msh"); - smm.applyBC(BC::Dirichlet::FixedValue(0., _y), "bottom"); - smm.applyBC(BC::Dirichlet::FixedValue(0., _x), "left"); + SolidMechanicsModel model(mesh); + model.initFull(_analysis_method = _static); - smm.setBaseName( "square_notch"); - smm.addDumpFieldVector( "displacement"); - smm.addDumpFieldVector( "internal_force"); - smm.addDumpField( "stress"); - smm.addDumpField( "grad_u"); - smm.addDumpField( "damage"); - smm.addDumpField( "blocked_dofs"); - smm.dump(); + model.setBaseName("solid_one_element"); + model.addDumpField("stress"); + model.addDumpField("grad_u"); + model.addDumpField("displacement"); + model.dump(); - auto & smm_solver = smm.getNonLinearSolver(); - smm_solver.set("max_iterations", 1000); - smm_solver.set("threshold", 1e-8); - smm_solver.set("convergence_type", _scc_residual); - - // coupling of models - SolidPhaseCoupler coupler(smm, pfm); + UInt nbSteps = 1000; - UInt nbSteps = 1500; - Real increment = 1.e-5; - - for (UInt s = 1; s < nbSteps; ++s) { - smm.applyBC(BC::Dirichlet::IncrementValue(increment, _y), "top"); - //smm.solveStep(); - coupler.solve(); - smm.dump(); - std::cout << "Step " << s << "/" << nbSteps << std::endl; - + for (UInt s = 0; s < nbSteps; ++s) { + applyDisplacement(model); + model.solveStep(); + model.dump(); } - + finalize(); + return EXIT_SUCCESS; + +} + + +/* -------------------------------------------------------------------------- */ +void applyDisplacement(SolidMechanicsModel & model) { + auto & displacement = model.getDisplacement(); + + auto & positions = model.getMesh().getNodes(); + auto & blocked_dofs = model.getBlockedDOFs(); + + + for (UInt n = 0; n < model.getMesh().getNbNodes(); ++n) { + if (positions(n, 1) == -0.5) { + displacement(n, 0) = 0; + displacement(n, 1) = 0; + blocked_dofs(n, 0) = true; + blocked_dofs(n ,1) = true; + } + else { + displacement(n, 0) = 0; + displacement(n, 1) += 1.e-4; + blocked_dofs(n, 0) = true; + blocked_dofs(n ,1) = true; + } + } } diff --git a/test/test_model/test_phase_field_model/test_square2d.geo b/test/test_model/test_phase_field_model/test_square2d.geo new file mode 100644 index 000000000..7d80bcfb0 --- /dev/null +++ b/test/test_model/test_phase_field_model/test_square2d.geo @@ -0,0 +1,29 @@ +// Mesh size +h = 0.05; + +// Dimensions of the square +Lx = 1.0; +Ly = 1.0; + +// ------------------------------------------ +// Geometry +// ------------------------------------------ +Point(1) = { 0.0, 0.0, 0.0, h}; +Point(2) = { Lx, 0.0, 0.0, h}; +Point(3) = { Lx, Ly, 0.0, h}; +Point(4) = { 0.0, Ly, 0.0, h}; + +Line(1) = {1, 2}; +Line(2) = {2, 3}; +Line(3) = {3, 4}; +Line(4) = {4, 1}; + +Line Loop(1) = {1:4}; + +Plane Surface(1) = {1}; + +Physical Surface(1) = {1}; +Physical Line("Crack") = {1}; +Physical Line("Fixed_x") = {4}; +Physical Line("Traction") = {2}; +Physical Line("Free") = {3};