diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc
index fd4533968..c50187f91 100644
--- a/src/model/contact_mechanics/contact_mechanics_model.cc
+++ b/src/model/contact_mechanics/contact_mechanics_model.cc
@@ -1,591 +1,596 @@
/**
* @file coontact_mechanics_model.cc
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Tue May 08 2012
* @date last modification: Wed Feb 21 2018
*
* @brief Contact mechanics 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "contact_mechanics_model.hh"
#include "boundary_condition_functor.hh"
#include "dumpable_inline_impl.hh"
#include "integrator_gauss.hh"
#include "shape_lagrange.hh"
#ifdef AKANTU_USE_IOHELPER
#include "dumper_iohelper_paraview.hh"
#endif
/* -------------------------------------------------------------------------- */
#include <algorithm>
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
ContactMechanicsModel::ContactMechanicsModel(
Mesh & mesh, Array<Real> & positions, UInt dim, const ID & id,
const MemoryID & memory_id, std::shared_ptr<DOFManager> dof_manager,
const ModelType model_type)
: Model(mesh, model_type, dof_manager, dim, id, memory_id),
current_positions(positions) {
AKANTU_DEBUG_IN();
this->registerFEEngineObject<MyFEEngineType>("ContactMechanicsModel", mesh,
Model::spatial_dimension);
#if defined(AKANTU_USE_IOHELPER)
this->mesh.registerDumper<DumperParaview>("contact_mechanics", id, true);
this->mesh.addDumpMeshToDumper("contact_mechanics", mesh,
Model::spatial_dimension, _not_ghost,
_ek_regular);
#endif
this->registerDataAccessor(*this);
this->detector = std::make_unique<ContactDetector>(this->mesh, current_positions,
id + ":contact_detector");
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
ContactMechanicsModel::ContactMechanicsModel(
Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id,
std::shared_ptr<DOFManager> dof_manager, const ModelType model_type)
: ContactMechanicsModel(mesh, mesh.getNodes(), dim, id, memory_id,
dof_manager, model_type) {}
/* -------------------------------------------------------------------------- */
ContactMechanicsModel::~ContactMechanicsModel() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::initFullImpl(const ModelOptions & options) {
Model::initFullImpl(options);
// initalize the resolutions
if (this->parser.getLastParsedFile() != "") {
this->instantiateResolutions();
}
this->initResolutions();
this->initBC(*this, *displacement, *displacement_increment, *external_force);
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::instantiateResolutions() {
ParserSection model_section;
bool is_empty;
std::tie(model_section, is_empty) = this->getParserSection();
if (not is_empty) {
auto model_resolutions =
model_section.getSubSections(ParserType::_contact_resolution);
for (const auto & section : model_resolutions) {
this->registerNewResolution(section);
}
}
auto sub_sections =
this->parser.getSubSections(ParserType::_contact_resolution);
for (const auto & section : sub_sections) {
this->registerNewResolution(section);
}
if (resolutions.empty())
AKANTU_EXCEPTION("No contact resolutions where instantiated for the model"
<< getID());
are_resolutions_instantiated = true;
}
/* -------------------------------------------------------------------------- */
Resolution &
ContactMechanicsModel::registerNewResolution(const ParserSection & section) {
std::string res_name;
std::string res_type = section.getName();
std::string opt_param = section.getOption();
try {
std::string tmp = section.getParameter("name");
res_name = tmp; /** this can seem weird, but there is an ambiguous operator
* overload that i couldn't solve. @todo remove the
* weirdness of this code
*/
} catch (debug::Exception &) {
AKANTU_ERROR("A contact resolution of type \'"
<< res_type
<< "\' in the input file has been defined without a name!");
}
Resolution & res = this->registerNewResolution(res_name, res_type, opt_param);
res.parseSection(section);
return res;
}
/* -------------------------------------------------------------------------- */
Resolution & ContactMechanicsModel::registerNewResolution(
const ID & res_name, const ID & res_type, const ID & opt_param) {
AKANTU_DEBUG_ASSERT(resolutions_names_to_id.find(res_name) ==
resolutions_names_to_id.end(),
"A resolution with this name '"
<< res_name << "' has already been registered. "
<< "Please use unique names for resolutions");
UInt res_count = resolutions.size();
resolutions_names_to_id[res_name] = res_count;
std::stringstream sstr_res;
sstr_res << this->id << ":" << res_count << ":" << res_type;
ID res_id = sstr_res.str();
std::unique_ptr<Resolution> resolution =
ResolutionFactory::getInstance().allocate(res_type, spatial_dimension,
opt_param, *this, res_id);
resolutions.push_back(std::move(resolution));
return *(resolutions.back());
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::initResolutions() {
AKANTU_DEBUG_ASSERT(resolutions.size() != 0,
"No resolutions to initialize !");
if (!are_resolutions_instantiated)
instantiateResolutions();
// \TODO check if each resolution needs a initResolution() method
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::initModel() {
getFEEngine().initShapeFunctions(_not_ghost);
getFEEngine().initShapeFunctions(_ghost);
}
/* -------------------------------------------------------------------------- */
FEEngine & ContactMechanicsModel::getFEEngineBoundary(const ID & name) {
return dynamic_cast<FEEngine &>(
getFEEngineClassBoundary<MyFEEngineType>(name));
}
/* -------------------------------------------------------------------------- */
-void ContactMechanicsModel::initSolver(
- TimeStepSolverType /*time_step_solver_type*/, NonLinearSolverType) {
+void ContactMechanicsModel::initSolver(TimeStepSolverType time_step_solver_type,
+ NonLinearSolverType) {
+
+ auto & dof_manager = this->getDOFManager();
+
+ /* -------------------------------------------------------------------------- */
+ // for alloc type of solvers
this->allocNodalField(this->displacement, spatial_dimension, "displacement");
this->allocNodalField(this->displacement_increment, spatial_dimension,
"displacement_increment");
- this->allocNodalField(this->contact_force, spatial_dimension,
- "contact_force");
+ this->allocNodalField(this->internal_force, spatial_dimension,
+ "internal_force");
this->allocNodalField(this->external_force, spatial_dimension,
"external_force");
this->allocNodalField(this->normals, spatial_dimension, "normals");
this->allocNodalField(this->gaps, 1, "gaps");
this->allocNodalField(this->nodal_area, 1, "areas");
this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs");
}
/* -------------------------------------------------------------------------- */
std::tuple<ID, TimeStepSolverType>
ContactMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) {
switch (method) {
case _explicit_contact: {
return std::make_tuple("explicit_contact", _tsst_dynamic);
}
case _implicit_contact: {
return std::make_tuple("implicit_contact", _tsst_static);
}
default:
return std::make_tuple("unkown", _tsst_not_defined);
}
}
/* -------------------------------------------------------------------------- */
ModelSolverOptions ContactMechanicsModel::getDefaultSolverOptions(
const TimeStepSolverType & type) const {
ModelSolverOptions options;
switch (type) {
case _tsst_dynamic: {
options.non_linear_solver_type = _nls_newton_raphson;
options.integration_scheme_type["displacement"] = _ist_pseudo_time;
options.solution_type["displacement"] = IntegrationScheme::_not_defined;
break;
}
case _tsst_static: {
options.non_linear_solver_type = _nls_newton_raphson;
options.integration_scheme_type["displacement"] = _ist_pseudo_time;
options.solution_type["displacement"] = IntegrationScheme::_not_defined;
break;
}
default:
AKANTU_EXCEPTION(type << " is not a valid time step solver type");
break;
}
return options;
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::assembleResidual() {
AKANTU_DEBUG_IN();
/* ------------------------------------------------------------------------ */
// computes the internal forces
this->assembleInternalForces();
+ /* ------------------------------------------------------------------------ */
+ this->getDOFManager().assembleToResidual("displacement", *this->internal_force,
+ 1);
+
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::assembleInternalForces() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the contact forces");
// assemble the forces due to local stresses
auto assemble = [&](auto && ghost_type) {
for (auto & resolution : resolutions) {
resolution->assembleInternalForces(ghost_type);
}
};
AKANTU_DEBUG_INFO("Assemble residual for local elements");
assemble(_not_ghost);
// assemble the stresses due to ghost elements
AKANTU_DEBUG_INFO("Assemble residual for ghost elements");
assemble(_ghost);
- /* ------------------------------------------------------------------------ */
- this->getDOFManager().assembleToResidual("displacement", *this->contact_force,
- 1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::search() {
this->detector->search(this->contact_map);
this->assembleFieldsFromContactMap();
this->computeNodalAreas<Surface::slave>();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::search(Array<Real> & increment) {
this->detector->search(this->contact_map);
for (auto & entry : contact_map) {
auto & element = entry.second;
const auto & connectivity = element.connectivity;
auto master_node = connectivity[1];
Vector<Real> u(spatial_dimension);
for (UInt s : arange(spatial_dimension)) {
u(s) = increment(master_node, s);
}
u *= -1.0;
const auto & normal = element.normal;
Real uv = Math::vectorDot(u.storage(), normal.storage(), spatial_dimension);
if (uv + element.gap <= 0) {
element.gap = abs(uv + element.gap);
} else {
element.gap = 0.0;
}
}
this->assembleFieldsFromContactMap();
this->computeNodalAreas<Surface::slave>();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::assembleFieldsFromContactMap() {
if (this->contact_map.empty())
AKANTU_ERROR(
"Contact map is empty, Please run search before assembling the fields");
for (auto & entry : contact_map) {
const auto & element = entry.second;
auto connectivity = element.connectivity;
auto node = connectivity(0);
(*gaps)[node] = element.gap;
for (UInt i = 0; i < spatial_dimension; ++i)
(*normals)(node, i) = element.normal[i];
}
}
/* -------------------------------------------------------------------------- */
template <Surface id> void ContactMechanicsModel::computeNodalAreas() {
this->nodal_area->clear();
this->external_force->clear();
this->applyBC(
BC::Neumann::FromHigherDim(Matrix<Real>::eye(spatial_dimension, 1)),
this->detector->getSurfaceId<id>());
for (auto && tuple :
zip(*nodal_area, make_view(*external_force, spatial_dimension))) {
auto & area = std::get<0>(tuple);
auto & force = std::get<1>(tuple);
for (auto & f : force)
area += pow(f, 2);
area = sqrt(area);
}
this->external_force->clear();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::beforeSolveStep() {
this->search();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::afterSolveStep() {}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::printself(std::ostream & stream, int indent) const {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
stream << space << "Contact Mechanics 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 << " + resolutions [" << std::endl;
for (auto & resolution : resolutions) {
resolution->printself(stream, indent + 1);
}
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
MatrixType ContactMechanicsModel::getMatrixType(const ID & matrix_id) {
if (matrix_id == "K")
return _symmetric;
return _mt_not_defined;
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::assembleMatrix(const ID & matrix_id) {
if (matrix_id == "K") {
this->assembleStiffnessMatrix();
}
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::assembleStiffnessMatrix() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the new stiffness matrix");
if (!this->getDOFManager().hasMatrix("K")) {
this->getDOFManager().getNewMatrix("K", getMatrixType("K"));
}
- this->getDOFManager().clearMatrix("K");
+ //this->getDOFManager().clearMatrix("K");
for (auto & resolution : resolutions) {
resolution->assembleStiffnessMatrix(_not_ghost);
}
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::assembleLumpedMatrix(const ID & /*matrix_id*/) {
AKANTU_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
dumper::Field *
ContactMechanicsModel::createNodalFieldBool(const std::string & /*field_name*/,
const std::string & /*group_name*/,
bool /*padding_flag*/) {
dumper::Field * field = nullptr;
return field;
}
/* -------------------------------------------------------------------------- */
dumper::Field *
ContactMechanicsModel::createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
std::map<std::string, Array<Real> *> real_nodal_fields;
- real_nodal_fields["contact_force"] = this->contact_force;
- real_nodal_fields["external_force"] = this->external_force;
+ real_nodal_fields["contact_force"] = this->internal_force;
real_nodal_fields["blocked_dofs"] = this->blocked_dofs;
real_nodal_fields["normals"] = this->normals;
real_nodal_fields["gaps"] = this->gaps;
real_nodal_fields["areas"] = this->nodal_area;
dumper::Field * field = nullptr;
if (padding_flag)
field = this->mesh.createNodalField(real_nodal_fields[field_name],
group_name, 3);
else
field =
this->mesh.createNodalField(real_nodal_fields[field_name], group_name);
return field;
}
#else
/* -------------------------------------------------------------------------- */
dumper::Field * ContactMechanicsModel::createNodalFieldReal(
__attribute__((unused)) const std::string & field_name,
__attribute__((unused)) const std::string & group_name,
__attribute__((unused)) bool padding_flag) {
return nullptr;
}
#endif
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::dump(const std::string & dumper_name) {
mesh.dump(dumper_name);
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::dump(const std::string & dumper_name, UInt step) {
mesh.dump(dumper_name, step);
}
/* ------------------------------------------------------------------------- */
void ContactMechanicsModel::dump(const std::string & dumper_name, Real time,
UInt step) {
mesh.dump(dumper_name, time, step);
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::dump() { mesh.dump(); }
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::dump(UInt step) { mesh.dump(step); }
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::dump(Real time, UInt step) {
mesh.dump(time, step);
}
/* -------------------------------------------------------------------------- */
UInt ContactMechanicsModel::getNbData(
const Array<Element> & 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);
}
AKANTU_DEBUG_OUT();
return size;
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::packData(CommunicationBuffer & /*buffer*/,
const Array<Element> & /*elements*/,
const SynchronizationTag & /*tag*/) const {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::unpackData(CommunicationBuffer & /*buffer*/,
const Array<Element> & /*elements*/,
const SynchronizationTag & /*tag*/) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
UInt ContactMechanicsModel::getNbData(
const Array<UInt> & dofs, const SynchronizationTag & /*tag*/) const {
AKANTU_DEBUG_IN();
UInt size = 0;
// UInt nb_nodes = mesh.getNbNodes();
AKANTU_DEBUG_OUT();
return size * dofs.size();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::packData(CommunicationBuffer & /*buffer*/,
const Array<UInt> & /*dofs*/,
const SynchronizationTag & /*tag*/) const {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::unpackData(CommunicationBuffer & /*buffer*/,
const Array<UInt> & /*dofs*/,
const SynchronizationTag & /*tag*/) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
} // namespace akantu
diff --git a/src/model/contact_mechanics/contact_mechanics_model.hh b/src/model/contact_mechanics/contact_mechanics_model.hh
index 181d78cba..9f5b43973 100644
--- a/src/model/contact_mechanics/contact_mechanics_model.hh
+++ b/src/model/contact_mechanics/contact_mechanics_model.hh
@@ -1,325 +1,325 @@
/**
* @file contact_mechanics_model.hh
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Tue Jul 27 2010
* @date last modification: Wed Feb 21 2018
*
* @brief Model of Contact Mechanics
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "boundary_condition.hh"
#include "contact_detector.hh"
#include "data_accessor.hh"
#include "fe_engine.hh"
#include "model.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_CONTACT_MECHANICS_MODEL_HH__
#define __AKANTU_CONTACT_MECHANICS_MODEL_HH__
namespace akantu {
class Resolution;
template <ElementKind kind, class IntegrationOrderFunctor>
class IntegratorGauss;
template <ElementKind kind> class ShapeLagrange;
} // namespace akantu
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
class ContactMechanicsModel : public Model,
public DataAccessor<Element>,
public DataAccessor<UInt>,
public BoundaryCondition<ContactMechanicsModel> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
using MyFEEngineType = FEEngineTemplate<IntegratorGauss, ShapeLagrange>;
public:
ContactMechanicsModel(
Mesh & mesh, UInt spatial_dimension = _all_dimensions,
const ID & id = "contact_mechanics_model", const MemoryID & memory_id = 0,
std::shared_ptr<DOFManager> dof_manager = nullptr,
const ModelType model_type = ModelType::_contact_mechanics_model);
ContactMechanicsModel(
Mesh & mesh, Array<Real> & positions,
UInt spatial_dimension = _all_dimensions,
const ID & id = "contact_mechanics_model", const MemoryID & memory_id = 0,
std::shared_ptr<DOFManager> dof_manager = nullptr,
const ModelType model_type = ModelType::_contact_mechanics_model);
~ContactMechanicsModel() override;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
protected:
/// initialize completely the model
void initFullImpl(const ModelOptions & options) override;
/// allocate all vectors
void initSolver(TimeStepSolverType, NonLinearSolverType) override;
/// initialize all internal arrays for resolutions
void initResolutions();
/// initialize the modelType
void initModel() override;
/// call back for the solver, computes the force residual
void assembleResidual() override;
/// get the type of matrix needed
MatrixType getMatrixType(const ID & matrix_id) override;
/// callback for the solver, this assembles different matrices
void assembleMatrix(const ID & matrix_id) override;
/// callback for the solver, this assembles the stiffness matrix
void assembleLumpedMatrix(const ID & matrix_id) override;
/// get some default values for derived classes
std::tuple<ID, TimeStepSolverType>
getDefaultSolverID(const AnalysisMethod & method) override;
ModelSolverOptions
getDefaultSolverOptions(const TimeStepSolverType & type) const;
/// 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;
/// function to print the containt of the class
void printself(std::ostream & stream, int indent = 0) const override;
/* ------------------------------------------------------------------------ */
/* Contact Detection */
/* ------------------------------------------------------------------------ */
public:
void search();
void search(Array<Real> &);
template <Surface id> void computeNodalAreas();
void assembleFieldsFromContactMap();
/* ------------------------------------------------------------------------ */
/* Contact Resolution */
/* ------------------------------------------------------------------------ */
public:
/// register an empty contact resolution of a given type
Resolution & registerNewResolution(const ID & res_name, const ID & res_type,
const ID & opt_param);
protected:
/// register a resolution in the dynamic database
Resolution & registerNewResolution(const ParserSection & res_section);
/// read the resolution files to instantiate all the resolutions
void instantiateResolutions();
/* ------------------------------------------------------------------------ */
/* Solver Interface */
/* ------------------------------------------------------------------------ */
public:
/// assembles the contact stiffness matrix
virtual void assembleStiffnessMatrix();
/// assembles the contant internal forces
virtual void assembleInternalForces();
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
FEEngine & getFEEngineBoundary(const ID & name = "") override;
/* ------------------------------------------------------------------------ */
/* 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;
void dump() override;
virtual void dump(UInt step);
virtual void dump(Real time, UInt step);
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);
/* ------------------------------------------------------------------------ */
/* Data Accessor inherited members */
/* ------------------------------------------------------------------------ */
public:
UInt getNbData(const Array<Element> & elements,
const SynchronizationTag & tag) const override;
void packData(CommunicationBuffer & buffer, const Array<Element> & elements,
const SynchronizationTag & tag) const override;
void unpackData(CommunicationBuffer & buffer, const Array<Element> & elements,
const SynchronizationTag & tag) override;
UInt getNbData(const Array<UInt> & dofs,
const SynchronizationTag & tag) const override;
void packData(CommunicationBuffer & buffer, const Array<UInt> & dofs,
const SynchronizationTag & tag) const override;
void unpackData(CommunicationBuffer & buffer, const Array<UInt> & dofs,
const SynchronizationTag & tag) override;
protected:
/// contact detection class
friend class ContactDetector;
/// contact resolution class
friend class Resolution;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// return the dimension of the system space
AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt);
/// get the ContactMechanicsModel::displacement vector
AKANTU_GET_MACRO(Displacement, *displacement, Array<Real> &);
/// get the ContactMechanicsModel::increment vector \warn only consistent
/// if ContactMechanicsModel::setIncrementFlagOn has been called before
AKANTU_GET_MACRO(Increment, *displacement_increment, Array<Real> &);
- /// get the ContactMechanics::contact_force vector (internal forces)
- AKANTU_GET_MACRO(InternalForce, *contact_force, Array<Real> &);
+ /// get the ContactMechanics::internal_force vector (internal forces)
+ AKANTU_GET_MACRO(InternalForce, *internal_force, Array<Real> &);
/// get the ContactMechanicsModel::external_force vector (external forces)
AKANTU_GET_MACRO(ExternalForce, *external_force, Array<Real> &);
/// get the ContactMechanicsModel::force vector (external forces)
Array<Real> & getForce() {
AKANTU_DEBUG_WARNING("getForce was maintained for backward compatibility, "
"use getExternalForce instead");
return *external_force;
}
/// get the ContactMechanics::blocked_dofs vector
AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array<Real> &);
/// get the ContactMechanics::gaps (contact gaps)
AKANTU_GET_MACRO(Gaps, *gaps, Array<Real> &);
/// get the ContactMechanics::normals (normals on slave nodes)
AKANTU_GET_MACRO(Normals, *normals, Array<Real> &);
/// get the ContactMechanics::areas (nodal areas)
AKANTU_GET_MACRO(NodalArea, *nodal_area, Array<Real> &);
- /// get the ContactMechanics::contact_force vector (internal forces)
+ /// get the ContactMechanics::internal_force vector (internal forces)
AKANTU_GET_MACRO(CurrentPositions, current_positions, Array<Real> &);
/// get the contat map
inline std::map<UInt, ContactElement> & getContactMap() {
return contact_map;
}
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// tells if the resolutions are instantiated
bool are_resolutions_instantiated;
/// displacements array
Array<Real> * displacement;
/// increment of displacement
Array<Real> * displacement_increment{nullptr};
/// contact forces array
- Array<Real> * contact_force{nullptr};
+ Array<Real> * internal_force{nullptr};
/// external forces array
Array<Real> * external_force{nullptr};
/// boundary vector
Array<Real> * blocked_dofs{nullptr};
/// gaps
Array<Real> * gaps{nullptr};
/// normals
Array<Real> * normals{nullptr};
/// tangents
Array<Real> * tangents{nullptr};
/// nodal areas
Array<Real> * nodal_area{nullptr};
/// array of current position used during update residual
Array<Real> & current_positions;
/// contact detection
std::unique_ptr<ContactDetector> detector;
/// list of contact resolutions
std::vector<std::unique_ptr<Resolution>> resolutions;
/// mapping between resolution name and resolution internal id
std::map<std::string, UInt> resolutions_names_to_id;
/// mapping between slave node its respective contact element
std::map<UInt, ContactElement> contact_map;
};
} // namespace akantu
/* ------------------------------------------------------------------------ */
/* inline functions */
/* ------------------------------------------------------------------------ */
#include "parser.hh"
#include "resolution.hh"
/* ------------------------------------------------------------------------ */
#endif /* __AKANTU_CONTACT_MECHANICS_MODEL_HH__ */
diff --git a/src/model/contact_mechanics/resolution.cc b/src/model/contact_mechanics/resolution.cc
index 85e718cc3..c3ce6db6f 100644
--- a/src/model/contact_mechanics/resolution.cc
+++ b/src/model/contact_mechanics/resolution.cc
@@ -1,356 +1,357 @@
/**
* @file resolution.cc
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Mon Jan 7 2019
* @date last modification: Mon Jan 7 2019
*
* @brief Implementation of common part of the contact resolution 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "resolution.hh"
#include "contact_mechanics_model.hh"
#include "sparse_matrix.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
Resolution::Resolution(ContactMechanicsModel & model, const ID & id)
: Memory(id, model.getMemoryID()), Parsable(ParserType::_contact_resolution, id),
fem(model.getFEEngine()),
name(""), model(model),
spatial_dimension(model.getMesh().getSpatialDimension()){
AKANTU_DEBUG_IN();
this->initialize();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Resolution::~Resolution() = default;
/* -------------------------------------------------------------------------- */
void Resolution::initialize() {
registerParam("name", name, std::string(), _pat_parsable | _pat_readable);
registerParam("mu", mu, Real(0.), _pat_parsable | _pat_modifiable,
"Friciton Coefficient");
}
/* -------------------------------------------------------------------------- */
void Resolution::printself(std::ostream & stream, int indent) const {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
std::string type = getID().substr(getID().find_last_of(':') + 1);
stream << space << "Contact Resolution " << type << " [" << std::endl;
Parsable::printself(stream, indent);
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
void Resolution::assembleInternalForces(GhostType /*ghost_type*/) {
AKANTU_DEBUG_IN();
auto & internal_force =
const_cast<Array<Real> &>(model.getInternalForce());
const auto local_nodes =
model.getMesh().getElementGroup(name).getNodes();
auto & nodal_area =
const_cast<Array<Real> &>(model.getNodalArea());
auto & contact_map = model.getContactMap();
for (auto & slave: local_nodes) {
if (contact_map.find(slave) == contact_map.end())
continue;
auto & element = contact_map[slave];
const auto & conn = element.connectivity;
const auto & type = element.master.type;
auto nb_nodes_master = Mesh::getNbNodesPerElement(type);
Vector<Real> shapes(nb_nodes_master);
Matrix<Real> dnds(spatial_dimension - 1, nb_nodes_master);
#define GET_SHAPES_NATURAL(type) \
ElementClass<type>::computeShapes(element.projection, shapes)
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPES_NATURAL);
#undef GET_SHAPES_NATURAL
#define GET_SHAPE_DERIVATIVES_NATURAL(type) \
ElementClass<type>::computeDNDS(element.projection, dnds)
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL);
#undef GET_SHAPE_DERIVATIVES_NATURAL
Vector<Real> fc(conn.size() * spatial_dimension);
- Matrix<Real> surface_matrix(spatial_dimension - 1, spatial_dimension - 1);
- computeSurfaceMatrix(element.tangents, surface_matrix);
+ Matrix<Real> m_alpha_beta(spatial_dimension - 1, spatial_dimension - 1);
+ computeMetricTensor(element.tangents, m_alpha_beta);
Vector<Real> n(conn.size() * spatial_dimension);
computeN(n, shapes, element.normal);
computeNormalForce(fc, n, element.gap);
Array<Real> t_alpha(conn.size() * spatial_dimension, spatial_dimension - 1);
Array<Real> n_alpha(conn.size() * spatial_dimension, spatial_dimension - 1);
Array<Real> d_alpha(conn.size() * spatial_dimension, spatial_dimension - 1);
computeTalpha(t_alpha, shapes, element.tangents);
computeNalpha(n_alpha, dnds, element.normal);
- computeDalpha(d_alpha, n_alpha, t_alpha, surface_matrix, element.gap);
+ computeDalpha(d_alpha, n_alpha, t_alpha, m_alpha_beta, element.gap);
//computeFrictionForce(fc, d_alpha, gap);
UInt nb_degree_of_freedom = internal_force.getNbComponent();
for (UInt i = 0; i < conn.size(); ++i) {
UInt n = conn[i];
for (UInt j = 0; j < nb_degree_of_freedom; ++j) {
UInt offset_node = n * nb_degree_of_freedom + j;
internal_force[offset_node] += fc[i*nb_degree_of_freedom + j];
internal_force[offset_node] *= nodal_area[n];
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Resolution::assembleStiffnessMatrix(GhostType /*ghost_type*/) {
AKANTU_DEBUG_IN();
auto & stiffness =
const_cast<SparseMatrix &>(model.getDOFManager().getMatrix("K"));
const auto local_nodes =
model.getMesh().getElementGroup(name).getNodes();
auto & nodal_area =
const_cast<Array<Real> &>(model.getNodalArea());
auto & contact_map = model.getContactMap();
for (auto & slave: local_nodes) {
if (contact_map.find(slave) == contact_map.end()) {
continue;
}
auto & master = contact_map[slave].master;
auto & gap = contact_map[slave].gap;
auto & projection = contact_map[slave].projection;
auto & normal = contact_map[slave].normal;
const auto & connectivity = contact_map[slave].connectivity;
const ElementType & type = master.type;
UInt nb_nodes_master = Mesh::getNbNodesPerElement(master.type);
Vector<Real> shapes(nb_nodes_master);
Matrix<Real> shapes_derivatives(spatial_dimension - 1, nb_nodes_master);
#define GET_SHAPES_NATURAL(type) \
ElementClass<type>::computeShapes(projection, shapes)
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPES_NATURAL);
#undef GET_SHAPES_NATURAL
#define GET_SHAPE_DERIVATIVES_NATURAL(type) \
ElementClass<type>::computeDNDS(projection, shapes_derivatives)
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL);
#undef GET_SHAPE_DERIVATIVES_NATURAL
Matrix<Real> elementary_stiffness(connectivity.size() * spatial_dimension,
connectivity.size() * spatial_dimension);
Matrix<Real> tangents(spatial_dimension - 1, spatial_dimension);
Matrix<Real> global_coords(nb_nodes_master, spatial_dimension);
computeCoordinates(master, global_coords);
computeTangents(shapes_derivatives, global_coords, tangents);
- Matrix<Real> surface_matrix(spatial_dimension - 1, spatial_dimension - 1);
- computeSurfaceMatrix(tangents, surface_matrix);
+ Matrix<Real> m_alpha_beta(spatial_dimension - 1, spatial_dimension - 1);
+ computeMetricTensor(tangents, m_alpha_beta);
Vector<Real> n(connectivity.size() * spatial_dimension);
Array<Real> t_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1);
Array<Real> n_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1);
Array<Real> d_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1);
computeN( n, shapes, normal);
computeTalpha( t_alpha, shapes, tangents);
computeNalpha( n_alpha, shapes_derivatives, normal);
- computeDalpha( d_alpha, n_alpha, t_alpha, surface_matrix, gap);
-
+ computeDalpha( d_alpha, n_alpha, t_alpha, m_alpha_beta, gap);
+
+ /*
Array<Real> t_alpha_beta(conn.size() * spatial_dimension, (spatial_dimension - 1) * (spatial_dimension -1));
Array<Real> n_alpha_beta(conn.size() * spatial_dimension, (spatial_dimension - 1) * (spatial_dimension -1));
Array<Real> p_alpha(conn.size() * spatial_dimension, spatial_dimension - 1);
-
+ */
Matrix<Real> kc(connectivity.size() * spatial_dimension,
connectivity.size() * spatial_dimension);
- computeTangentModuli(kc, n, n_alpha, d_alpha, surface_matrix, gap);
+ computeTangentModuli(kc, n, n_alpha, d_alpha, m_alpha_beta, gap);
std::vector<UInt> equations;
UInt nb_degree_of_freedom = model.getSpatialDimension();
std::vector<Real> areas;
for (UInt i : arange(connectivity.size())) {
UInt n = connectivity[i];
for (UInt j : arange(nb_degree_of_freedom)) {
equations.push_back(n * nb_degree_of_freedom + j);
areas.push_back(nodal_area[n]);
}
}
for (UInt i : arange(kc.rows())) {
UInt row = equations[i];
for (UInt j : arange(kc.cols())) {
UInt col = equations[j];
kc(i, j) *= areas[i];
stiffness.add(row, col, kc(i, j));
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Resolution::computeTangents(Matrix<Real> & shapes_derivatives, Matrix<Real> & global_coords,
Matrix<Real> & tangents) {
tangents.mul<false, false>(shapes_derivatives, global_coords);
}
/* -------------------------------------------------------------------------- */
-void Resolution::computeSurfaceMatrix(Matrix<Real> & tangents, Matrix<Real> & surface_matrix) {
+void Resolution::computeMetricTensor(Matrix<Real> & tangents, Matrix<Real> & m_alpha_beta) {
- surface_matrix.mul<false, true>(tangents, tangents);
- surface_matrix = surface_matrix.inverse();
+ m_alpha_beta.mul<false, true>(tangents, tangents);
+ m_alpha_beta = m_alpha_beta.inverse();
}
/* -------------------------------------------------------------------------- */
void Resolution::computeN(Vector<Real> & n, Vector<Real> & shapes, Vector<Real> & normal) {
UInt dim = normal.size();
for (UInt i = 0; i < dim; ++i) {
n[i] = normal[i];
for (UInt j = 0; j < shapes.size(); ++j) {
n[(1 + j) * dim + i] = -normal[i] * shapes[j];
}
}
}
/* -------------------------------------------------------------------------- */
void Resolution::computeTalpha(Array<Real> & t_alpha, Vector<Real> & shapes,
Matrix<Real> & tangents) {
t_alpha.clear();
for (auto && values:
zip(tangents.transpose(),
make_view(t_alpha, t_alpha.size()))) {
auto & tangent = std::get<0>(values);
auto & t_s = std::get<1>(values);
for (UInt i : arange(spatial_dimension)) {
t_s[i] = -tangent(i);
for (UInt j : arange(shapes.size())) {
t_s[(1 + j)*spatial_dimension + i] = -shapes[j] * tangent(i);
}
}
}
}
/* -------------------------------------------------------------------------- */
void Resolution::computeNalpha(Array<Real> & n_alpha, Matrix<Real> & shapes_derivatives,
Vector<Real> & normal) {
n_alpha.clear();
for (auto && values:
zip(shapes_derivatives.transpose(),
make_view(n_alpha, n_alpha.size()))) {
auto & dnds = std::get<0>(values);
auto & n_s = std::get<1>(values);
for (UInt i : arange(spatial_dimension)) {
n_s[i] = 0;
for (UInt j : arange(dnds.size())) {
n_s[(1 + j)*spatial_dimension + i] = -dnds(j)*normal[i];
}
}
}
}
/* -------------------------------------------------------------------------- */
void Resolution::computeDalpha(Array<Real> & d_alpha, Array<Real> & n_alpha,
- Array<Real> & t_alpha, Matrix<Real> & surface_matrix,
+ Array<Real> & t_alpha, Matrix<Real> & m_alpha_beta,
Real & gap) {
d_alpha.clear();
- for (auto && entry : zip(surface_matrix.transpose(),
+ for (auto && entry : zip(m_alpha_beta.transpose(),
make_view(d_alpha, d_alpha.size()))) {
auto & a_s = std::get<0>(entry);
auto & d_s = std::get<1>(entry);
for (auto && values :
zip(arange(t_alpha.size()),
make_view(t_alpha, t_alpha.size()),
make_view(n_alpha, n_alpha.size()))) {
auto & index = std::get<0>(values);
auto & t_s = std::get<1>(values);
auto & n_s = std::get<2>(values);
d_s += (t_s + gap * n_s);
d_s *= a_s(index);
}
}
}
/* -------------------------------------------------------------------------- */
void Resolution::computeCoordinates(const Element & el, Matrix<Real> & coords) {
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type);
Vector<UInt> connect = model.getMesh().getConnectivity(el.type, _not_ghost)
.begin(nb_nodes_per_element)[el.element];
// change this to current position
auto & positions = model.getMesh().getNodes();
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
UInt node = connect[n];
for (UInt s: arange(spatial_dimension)) {
coords(n, s) = positions(node, s);
}
}
}
} // akantu
diff --git a/src/model/contact_mechanics/resolution.hh b/src/model/contact_mechanics/resolution.hh
index 8051d9428..67984d34a 100644
--- a/src/model/contact_mechanics/resolution.hh
+++ b/src/model/contact_mechanics/resolution.hh
@@ -1,230 +1,231 @@
/**
* @file resolution.hh
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Mon Jan 7 2019
* @date last modification: Mon Jan 7 2019
*
* @brief Mother class for all contact resolutions
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_factory.hh"
#include "aka_memory.hh"
#include "parsable.hh"
#include "parser.hh"
#include "fe_engine.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_RESOLUTION_HH__
#define __AKANTU_RESOLUTION_HH__
/* -------------------------------------------------------------------------- */
namespace akantu {
class Model;
class ContactMechanicsModel;
} // namespace akantu
namespace akantu {
/**
* Interface of all contact resolutions
* Prerequisites for a new resolution
* - inherit from this class
* - implement the following methods:
* \code
*
* virtual void computeNormalForce();
* virtual void computeFrictionForce();
*
* \endcode
*
*/
class Resolution : public Memory,
public Parsable {
/* ------------------------------------------------------------------------ */
/* Constructor/Destructor */
/* ------------------------------------------------------------------------ */
public:
/// instantiate contact resolution with defaults
Resolution(ContactMechanicsModel & model, const ID & id = "");
/// Destructor
~Resolution() override;
protected:
void initialize();
/// computes coordinates of a given element
void computeCoordinates(const Element & , Matrix<Real> &);
/// computes tangents
- void computeTangents(Matrix<Real> & /* shapes_derivatives */,
- Matrix<Real> & /* global_coords */,
- Matrix<Real> & /* tangents */);
+ void computeTangents(Matrix<Real> & ,
+ Matrix<Real> & ,
+ Matrix<Real> & );
- /// computes surface metric matrix
- void computeSurfaceMatrix(Matrix<Real> & /* tangents */,
- Matrix<Real> & /* surface_matrix */);
+ /// computes metric tensor (@f$m_{\alpha\beta}@f$) where @f$\alpha,
+ /// \beta@f$ are surface directions
+ void computeMetricTensor(Matrix<Real> &,
+ Matrix<Real> &);
/// computes N array
- void computeN(Vector<Real> & /* n */,
- Vector<Real> & /* shapes */,
- Vector<Real> & /* normal */);
+ void computeN(Vector<Real> &,
+ Vector<Real> &,
+ Vector<Real> &);
- /// computes N_{\alpha} where \alpha is number of surface dimensions
- void computeNalpha(Array<Real> & /* n_alpha */,
- Matrix<Real> & /* shapes_derivatives */,
- Vector<Real> & /* normal */);
+ /// computes (@f$N_{\alpha}@f$) where \alpha is surface dimension
+ void computeNalpha(Array<Real> &,
+ Matrix<Real> &,
+ Vector<Real> &);
- /// computes T_{\alpha} where \alpha is surface dimenion
+ /// computes (@f$T_{\alpha}@f$) where @f$\alpha@f$ is surface dimension
void computeTalpha(Array<Real> & ,
Vector<Real> & ,
Matrix<Real> & );
- /// computes D_{\alpha} where \alpha is number of surface dimensions
- void computeDalpha(Array<Real> & /* d_alpha */,
- Array<Real> & /* n_alpha */,
- Array<Real> & /* t_alpha */,
- Matrix<Real> & /* surface_matrix */,
- Real & /* gap */);
+ /// computes (@f$D_{\alpha}@f$) where @f$\alpha@f$ is surface dimension
+ void computeDalpha(Array<Real> & ,
+ Array<Real> & ,
+ Array<Real> & ,
+ Matrix<Real> &,
+ Real & );
/* ------------------------------------------------------------------------ */
/* Functions that resolutions can/should reimplement */
/* ------------------------------------------------------------------------ */
protected:
/// computes the normal force
virtual void computeNormalForce(__attribute__((unused)) Vector<Real> & /* force */,
__attribute__((unused)) Vector<Real> & /* n */,
__attribute__((unused)) Real & /* gap */) {
AKANTU_TO_IMPLEMENT();
}
/// computes the friction force
virtual void computeFrictionForce(__attribute__((unused)) Vector<Real> & /* force */,
__attribute__((unused)) Array<Real> & /* d_alpha */,
__attribute__((unused)) Real & /* gap */) {
AKANTU_TO_IMPLEMENT();
}
/// compute the tangent moduli
virtual void computeTangentModuli(__attribute__((unused)) Matrix<Real> &,
__attribute__((unused)) Vector<Real> & /* n */,
__attribute__((unused)) Array<Real> & /* n_alpha */,
__attribute__((unused)) Array<Real> & /* d_alpha */,
__attribute__((unused)) Matrix<Real> & /* A */,
__attribute__((unused)) Real & /* gap */) {
AKANTU_TO_IMPLEMENT();
}
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// assemble the residual for this resolution
void assembleInternalForces(GhostType ghost_type);
/// assemble the stiffness matrix for this resolution
void assembleStiffnessMatrix(GhostType ghost_type);
public:
/// function to print the contain of the class
void printself(std::ostream & stream, int indent = 0) const override;
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// Link to the fem object in the model
FEEngine & fem;
/// resolution name
std::string name;
/// model to which the resolution belong
ContactMechanicsModel & model;
/// friciton coefficient : mu
Real mu;
/// spatial dimension
UInt spatial_dimension;
};
/// standard output stream operator
inline std::ostream & operator<<(std::ostream & stream,
const Resolution & _this) {
_this.printself(stream);
return stream;
}
} // namespace akantu
/* -------------------------------------------------------------------------- */
namespace akantu {
using ResolutionFactory =
Factory<Resolution, ID, UInt, const ID &, ContactMechanicsModel &, const ID &>;
/// macaulay bracket to convert positive gap to zero
template <typename T>
T macaulay(T var) {return var < 0 ? 0 : var; }
template <typename T>
T heaviside(T var) {return var < 0 ? 0 : 1; }
} // namespace akantu
#define INSTANTIATE_RESOLUTION_ONLY(res_name) \
class res_name
#define RESOLUTION_DEFAULT_PER_DIM_ALLOCATOR(id, res_name) \
[](UInt dim, const ID &, ContactMechanicsModel & model, \
const ID & id) -> std::unique_ptr<Resolution> { \
switch (dim) { \
case 1: \
return std::make_unique<res_name>(model, id); \
case 2: \
return std::make_unique<res_name>(model, id); \
case 3: \
return std::make_unique<res_name>(model, id); \
default: \
AKANTU_EXCEPTION("The dimension " \
<< dim << "is not a valid dimension for the contact resolution " \
<< #id); \
} \
}
#define INSTANTIATE_RESOLUTION(id, res_name) \
INSTANTIATE_RESOLUTION_ONLY(res_name); \
static bool resolution_is_alocated_##id[[gnu::unused]] = \
ResolutionFactory::getInstance().registerAllocator( \
#id, RESOLUTION_DEFAULT_PER_DIM_ALLOCATOR(id, res_name))
#endif /* __AKANTU_RESOLUTION_HH__ */
diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty.cc b/src/model/contact_mechanics/resolutions/resolution_penalty.cc
index 18e3d19a2..eec5c5e12 100644
--- a/src/model/contact_mechanics/resolutions/resolution_penalty.cc
+++ b/src/model/contact_mechanics/resolutions/resolution_penalty.cc
@@ -1,154 +1,180 @@
/**
* @file resolution_penalty.cc
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Mon Jan 7 2019
* @date last modification: Mon Jan 7 2019
*
* @brief Specialization of the resolution class for the penalty method
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "resolution_penalty.hh"
namespace akantu {
/* -------------------------------------------------------------------------- */
ResolutionPenalty::ResolutionPenalty(ContactMechanicsModel & model,
const ID & id)
: Resolution(model, id) {
AKANTU_DEBUG_IN();
this->initialize();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::initialize() {
this->registerParam("epsilon", epsilon, Real(0.), _pat_parsable | _pat_modifiable,
"Normal penalty parameter");
this->registerParam("epsilon_t", epsilon_t, Real(0.), _pat_parsable | _pat_modifiable,
"Tangential penalty parameter");
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeNormalForce(Vector<Real> & force, Vector<Real> & n,
Real & gap) {
force.clear();
Real tn = gap * epsilon;
tn = macaulay(tn);
for (UInt i : arange(force.size())) {
force[i] += tn * n[i];
}
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeFrictionForce(Vector<Real> & force, Array<Real> & d_alpha,
Real & gap) {
Vector<Real> tractions(d_alpha.getNbComponent());
- computeFrictionalTraction(tractions);
+ /*computeFrictionalTraction(tractions);
for (auto && values:
zip(tractions,
make_view(d_alpha, d_alpha.size()))) {
auto & t_s = std::get<0>(values);
auto & d_s = std::get<1>(values);
force += d_s * t_s;
- }
+ }*/
}
+/* -------------------------------------------------------------------------- */
+void ResolutionPenalty::computeFrictionTraction(Real & gap) {
+
+ Real tn = gap * epsilon;
+ tn = macaulay(tn);
+
+ // delta_xi = xi- previous_xi
+ //trial_traction = previous_traction + epsilon * metric_tensor * delta_xi;
+ //trial_slip_function = trial_traction.mod - mu * tn;
+
+ // Stick condition
+ // if trial_slip_function <= 0
+ // traction = trial_traction;
+
+ // Slip condition
+ // else
+ // traction = mu * tn * trial_traction / trial_traction.norm()'
+
+
+}
+
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeTangentModuli(Matrix<Real> & kc, Vector<Real> & n,
Array<Real> & n_alpha, Array<Real> & d_alpha,
Matrix<Real> & surface_matrix,
Real & gap) {
computeNormalStiffness(kc, n, n_alpha, d_alpha, surface_matrix, gap);
computeFrictionalStiffness(n, n_alpha, d_alpha, gap);
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeNormalStiffness(Matrix<Real> & ke, Vector<Real> & n,
Array<Real> & n_alpha, Array<Real> & d_alpha,
Matrix<Real> & surface_matrix, Real & gap) {
Real tn = gap * epsilon;
tn = macaulay(tn);
Matrix<Real> n_mat(n.storage(), n.size(), 1);
ke.mul<false, true>(n_mat, n_mat);
ke *= epsilon * heaviside(gap);
for (auto && values:
zip(make_view(n_alpha, n_alpha.size()),
make_view(d_alpha, d_alpha.size()))) {
auto & n_s = std::get<0>(values);
auto & d_s = std::get<1>(values);
Matrix<Real> ns_mat(n_s.storage(), n_s.size(), 1);
Matrix<Real> ds_mat(d_s.storage(), d_s.size(), 1);
Matrix<Real> tmp1(n_s.size(), n_s.size());
tmp1.mul<false, true>(ns_mat, ds_mat);
Matrix<Real> tmp2(n_s.size(), n_s.size());
tmp1.mul<false, true>(ds_mat, ns_mat);
ke -= (tmp1 + tmp2) * tn;
}
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeFrictionalStiffness(Vector<Real> & n,
Array<Real> & n_alpha, Array<Real> & d_alpha,
Real & gap) {
-
+
+ computeCommonModuli();
+ computeStickModuli();
+ computeSlipModuli();
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeCommonModuli(Real & gap) {
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeStickModuli() {
+ // epsilon_t *
+
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeSlipModuli() {
}
INSTANTIATE_RESOLUTION(penalty, ResolutionPenalty);
} // akantu
diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty.hh b/src/model/contact_mechanics/resolutions/resolution_penalty.hh
index 21031c0f5..2650e63c0 100644
--- a/src/model/contact_mechanics/resolutions/resolution_penalty.hh
+++ b/src/model/contact_mechanics/resolutions/resolution_penalty.hh
@@ -1,113 +1,116 @@
/**
* @file contact_resolution_penalty.hh
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Mon Jan 29 2018
*
* @brief Material isotropic thermo-elastic
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "resolution.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_RESOLUTION_PENALTY_HH__
#define __AKANTU_RESOLUTION_PENALTY_HH__
namespace akantu {
class ResolutionPenalty : public Resolution {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
ResolutionPenalty(ContactMechanicsModel & model, const ID & id = "");
~ResolutionPenalty() override = default;
protected:
/// initialize the resolution
void initialize();
/// local computation of stifnness matrix due to normal stress
void computeNormalStiffness(Matrix<Real> & kc,
Vector<Real> & n,
Array<Real> & n_alpha,
Array<Real> & d_alpha,
Matrix<Real> & surface_matrix,
Real & gap);
/// local computation of stiffness matrix due to frictional stress
void computeFrictionalStiffness(Vector<Real> & n,
Array<Real> & n_alpha,
Array<Real> & d_alpha,
Real & gap);
/// local computation of direct stiffness matrix due to friction,
/// this matrix is common for both stick and slip part
void computeCommonModuli(Real & gap);
/// local computaion of stiffness matrix due to stick state
void computeStickModuli();
/// local computation of stiffness matrix due to slip state
void computeSlipModuli();
+
+ ///
+ void computeFrictionalTraction();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// local computation of tangent moduli
void computeTangentModuli(Matrix<Real> &,
Vector<Real> & /* N */,
Array<Real> & /* N_alpha */,
Array<Real> & /* D_alpha */,
Matrix<Real> & /* A matrix */,
Real & /* gap */
) override;
/// local computation of normal force
void computeNormalForce(Vector<Real> & /* force vector */,
Vector<Real> & /* n */,
Real & /* gap */) override;
/// local computation of friction force
void computeFrictionForce(Vector<Real> & /* force vector */,
Array<Real> & /* D_alpha */,
Real & /* gap */) override;
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// penalty parameter for normal stress
Real epsilon;
/// penalty parameter for tangential stress
Real epsilon_t;
};
} // akantu
#endif /* __AKANTU_RESOLUTION_PENALTY_HH__ */
diff --git a/src/model/model_couplers/coupler_solid_contact.cc b/src/model/model_couplers/coupler_solid_contact.cc
index 7aba79ba4..f6e63cbd2 100644
--- a/src/model/model_couplers/coupler_solid_contact.cc
+++ b/src/model/model_couplers/coupler_solid_contact.cc
@@ -1,282 +1,415 @@
/**
* @file coupler_solid_contact_explicit.cc
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Thu Jan 17 2019
* @date last modification: Thu Jan 17 2019
*
* @brief class for coupling of solid mechanics and conatct mechanics
* model in explicit
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "coupler_solid_contact.hh"
#include "dumpable_inline_impl.hh"
#include "integrator_gauss.hh"
#include "shape_lagrange.hh"
#ifdef AKANTU_USE_IOHELPER
#include "dumper_iohelper_paraview.hh"
#endif
/* -------------------------------------------------------------------------- */
namespace akantu {
-CouplerSolidContact::CouplerSolidContact(SolidMechanicsModel & solid,
- ContactMechanicsModel & contact,
- UInt dim, const ID & id,
+CouplerSolidContact::CouplerSolidContact(Mesh & mesh, UInt dim, const ID & id,
std::shared_ptr<DOFManager> dof_manager,
const ModelType model_type)
- : Model(solid.getMesh(), model_type, dof_manager, dim, id), solid(solid),
- contact(contact) {
+ : Model(mesh, model_type, dof_manager, dim, id) {
AKANTU_DEBUG_IN();
this->registerFEEngineObject<MyFEEngineType>(
- "CouplerSolidContact", solid.getMesh(), Model::spatial_dimension);
+ "CouplerSolidContact", mesh, Model::spatial_dimension);
#if defined(AKANTU_USE_IOHELPER)
this->mesh.registerDumper<DumperParaview>("coupler_solid_contact", id, true);
- this->mesh.addDumpMeshToDumper("coupler_solid_contact", solid.getMesh(),
+ this->mesh.addDumpMeshToDumper("coupler_solid_contact", mesh,
Model::spatial_dimension, _not_ghost,
_ek_regular);
#endif
this->registerDataAccessor(*this);
+
+ solid = new SolidMechanicsModel(mesh, Model::spatial_dimension,
+ "solid_mechanics_model",
+ 0,
+ this->dof_manager);
+ contact = new ContactMechanicsModel(mesh, Model::spatial_dimension,
+ "contact_mechanics_model",
+ 0,
+ this->dof_manager);
+
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
CouplerSolidContact::~CouplerSolidContact() {}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::initFullImpl(const ModelOptions & options) {
Model::initFullImpl(options);
this->initBC(*this, *displacement, *displacement_increment, *external_force);
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::initModel() {
+
getFEEngine().initShapeFunctions(_not_ghost);
getFEEngine().initShapeFunctions(_ghost);
+
}
/* -------------------------------------------------------------------------- */
FEEngine & CouplerSolidContact::getFEEngineBoundary(const ID & name) {
return dynamic_cast<FEEngine &>(
getFEEngineClassBoundary<MyFEEngineType>(name));
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::initSolver(TimeStepSolverType, NonLinearSolverType) {
DOFManager & dof_manager = this->getDOFManager();
- this->allocNodalField(this->displacement, spatial_dimension, "displacement");
+ /*this->allocNodalField(this->displacement, spatial_dimension, "displacement");
this->allocNodalField(this->displacement_increment, spatial_dimension,
- "displacement_increment");
+ "displacement_increment");
this->allocNodalField(this->external_force, spatial_dimension,
"external_force");
if (not dof_manager.hasDOFs("displacement")) {
dof_manager.registerDOFs("displacement", *this->displacement, _dst_nodal);
- }
+ }*/
}
/* -------------------------------------------------------------------------- */
std::tuple<ID, TimeStepSolverType>
CouplerSolidContact::getDefaultSolverID(const AnalysisMethod & method) {
switch (method) {
case _explicit_contact: {
return std::make_tuple("explicit_contact", _tsst_dynamic);
}
case _implicit_contact: {
return std::make_tuple("implicit_contact", _tsst_static);
}
default:
return std::make_tuple("unkown", _tsst_not_defined);
}
}
/* -------------------------------------------------------------------------- */
ModelSolverOptions CouplerSolidContact::getDefaultSolverOptions(
const TimeStepSolverType & type) const {
ModelSolverOptions options;
switch (type) {
case _tsst_dynamic: {
options.non_linear_solver_type = _nls_newton_raphson;
options.integration_scheme_type["displacement"] = _ist_pseudo_time;
options.solution_type["displacement"] = IntegrationScheme::_not_defined;
break;
}
case _tsst_static: {
options.non_linear_solver_type = _nls_newton_raphson;
options.integration_scheme_type["displacement"] = _ist_pseudo_time;
options.solution_type["displacement"] = IntegrationScheme::_not_defined;
break;
}
default:
AKANTU_EXCEPTION(type << " is not a valid time step solver type");
break;
}
return options;
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::assembleResidual() {
- solid.assembleInternalForces();
- contact.assembleInternalForces();
- this->coupleExternalForces();
+ solid->assembleInternalForces();
+ contact->assembleInternalForces();
+
+ auto & contact_force = contact->getInternalForce();
+ auto & external_force = solid->getExternalForce();
+ auto & internal_force = solid->getInternalForce();
+
+ /*auto & blocked_dofs = solid->getBlockedDOFs();
+
+ for (auto && values : zip(make_view(external_force),
+ make_view(contact_force),
+ make_view(blocked_dofs))) {
+ auto & f_ext = std::get<0>(values);
+ auto & f_con = std::get<1>(values);
+ auto & blocked = std::get<2>(values);
+
+ if (!blocked)
+ f_ext = f_con;
+ }*/
+
+
+ /* ------------------------------------------------------------------------ */
+ this->getDOFManager().assembleToResidual("displacement",
+ external_force, 1);
+ this->getDOFManager().assembleToResidual("displacement",
+ internal_force, 1);
+ this->getDOFManager().assembleToResidual("displacement",
+ contact_force, 1);
}
/* -------------------------------------------------------------------------- */
-void CouplerSolidContact::beforeSolveStep() { contact.search(); }
+void CouplerSolidContact::beforeSolveStep() {
+ contact->search();
+}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::afterSolveStep() {}
/* -------------------------------------------------------------------------- */
MatrixType CouplerSolidContact::getMatrixType(const ID & matrix_id) {
if (matrix_id == "K")
return _symmetric;
return _mt_not_defined;
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::assembleMatrix(const ID & matrix_id) {
- if (matrix_id == "K")
- contact.assembleStiffnessMatrix();
+ if (matrix_id == "K") {
+ solid->assembleStiffnessMatrix();
+ contact->assembleStiffnessMatrix();
+ }
- //auto & solid_tss = solid.getTimeStepSolver();
- // solid.assembleMatrix(matrix_id);
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::assembleLumpedMatrix(const ID & /*matrix_id*/) {
AKANTU_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::coupleExternalForces() {
- auto & contact_force = contact.getInternalForce();
- auto & external_force = solid.getExternalForce();
- auto & blocked_dofs = solid.getBlockedDOFs();
+ auto & contact_force = contact->getInternalForce();
+ auto & external_force = solid->getExternalForce();
+ auto & blocked_dofs = solid->getBlockedDOFs();
for (auto && values : zip(make_view(external_force), make_view(contact_force),
make_view(blocked_dofs))) {
auto & f_ext = std::get<0>(values);
auto & f_con = std::get<1>(values);
auto & blocked = std::get<2>(values);
if (!blocked)
f_ext = f_con;
}
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::coupleStiffnessMatrices() {
auto & contact_stiffness =
- const_cast<SparseMatrix &>(contact.getDOFManager().getMatrix("K"));
+ const_cast<SparseMatrix &>(contact->getDOFManager().getMatrix("K"));
auto & solid_stiffness =
- const_cast<SparseMatrix &>(solid.getDOFManager().getMatrix("K"));
+ const_cast<SparseMatrix &>(solid->getDOFManager().getMatrix("K"));
solid_stiffness.add(contact_stiffness);
}
/* -------------------------------------------------------------------------- */
// void CouplerSolidContact::solve() {
// search for contacts
// contact.search();
// can be handled by beforeSolveStep()
/// assemble contact forces
/// contact.assembleInternalForces();
/// assemble contact stiffness matrix
/// contact.assembleStiffnessMatrix();
/// can be handled by solveStep but no need to solve it
/// couple the external forces
// this->coupleExternalForces();
// assemble the internal forces solid mechanics model
// assemble the stiffness forces
// couple the contact mechanics model stiffness to solid mechanics
// this->coupleStiffnessMatrices();
// and solve the solid mechanics model with new stifffness anf
// residual
/// all the above steps for solid mehcanics should be handled by the
// solveStep method with solvercall back to see that stifffness
// matrices are coupled
// ContactSolver callback(solid, contact);
// contact.solveStep();
// solid.solveStep(callback);
//}
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
+/* -------------------------------------------------------------------------- */
+dumper::Field * CouplerSolidContact::createElementalField(
+ const std::string & field_name, const std::string & group_name,
+ bool padding_flag, const UInt & spatial_dimension,
+ const ElementKind & kind) {
+
+ dumper::Field * field = nullptr;
+
+ field = solid->createElementalField(field_name, group_name, padding_flag,
+ spatial_dimension, kind);
+
+ return field;
+}
+
+
+
/* -------------------------------------------------------------------------- */
dumper::Field *
CouplerSolidContact::createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
- std::map<std::string, Array<Real> *> real_nodal_fields;
- real_nodal_fields["displacement"] = this->displacement;
- real_nodal_fields["external_force"] = this->external_force;
-
dumper::Field * field = nullptr;
- if (padding_flag)
- field = this->mesh.createNodalField(real_nodal_fields[field_name],
- group_name, 3);
- else
- field =
- this->mesh.createNodalField(real_nodal_fields[field_name], group_name);
+
+ field = solid->createNodalFieldReal(field_name, group_name, padding_flag);
return field;
}
+/* -------------------------------------------------------------------------- */
+dumper::Field * CouplerSolidContact::createNodalFieldBool(
+ const std::string & field_name, const std::string & group_name,
+ __attribute__((unused)) bool padding_flag) {
+
+ dumper::Field * field = nullptr;
+ field = solid->createNodalFieldBool(field_name, group_name, padding_flag);
+ return field;
+}
+
#else
+
/* -------------------------------------------------------------------------- */
-dumper::Field * CouplerSolidContact::createNodalFieldReal(
- __attribute__((unused)) const std::string & field_name,
- __attribute__((unused)) const std::string & group_name,
- __attribute__((unused)) bool padding_flag) {
+dumper::Field * CouplerSolidContact::createElementalField(const std::string &,
+ const std::string &,
+ bool, const UInt &,
+ const ElementKind &) {
+ return nullptr;
+}
+/* ----------------------------------------------------------------------- */
+dumper::Field * CouplerSolidContact::createNodalFieldReal(const std::string &,
+ const std::string &,
+ bool) {
return nullptr;
}
+/*-------------------------------------------------------------------*/
+dumper::Field * CouplerSolidContact::createNodalFieldBool(const std::string &,
+ const std::string &,
+ bool) {
+ return nullptr;
+}
+
#endif
+/* -------------------------------------------------------------------------- */
+UInt CouplerSolidContact::getNbData(
+ const Array<Element> & 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);
+ }
+
+ AKANTU_DEBUG_OUT();
+ return size;
+}
+
+/* -------------------------------------------------------------------------- */
+void CouplerSolidContact::packData(CommunicationBuffer & /*buffer*/,
+ const Array<Element> & /*elements*/,
+ const SynchronizationTag & /*tag*/) const {
+ AKANTU_DEBUG_IN();
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+void CouplerSolidContact::unpackData(CommunicationBuffer & /*buffer*/,
+ const Array<Element> & /*elements*/,
+ const SynchronizationTag & /*tag*/) {
+ AKANTU_DEBUG_IN();
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+UInt CouplerSolidContact::getNbData(
+ const Array<UInt> & dofs, const SynchronizationTag & /*tag*/) const {
+ AKANTU_DEBUG_IN();
+
+ UInt size = 0;
+ // UInt nb_nodes = mesh.getNbNodes();
+
+ AKANTU_DEBUG_OUT();
+ return size * dofs.size();
+}
+
+/* -------------------------------------------------------------------------- */
+void CouplerSolidContact::packData(CommunicationBuffer & /*buffer*/,
+ const Array<UInt> & /*dofs*/,
+ const SynchronizationTag & /*tag*/) const {
+ AKANTU_DEBUG_IN();
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+void CouplerSolidContact::unpackData(CommunicationBuffer & /*buffer*/,
+ const Array<UInt> & /*dofs*/,
+ const SynchronizationTag & /*tag*/) {
+ AKANTU_DEBUG_IN();
+
+ AKANTU_DEBUG_OUT();
+}
+
+
} // namespace akantu
diff --git a/src/model/model_couplers/coupler_solid_contact.hh b/src/model/model_couplers/coupler_solid_contact.hh
index e2c13107e..fa16f370c 100644
--- a/src/model/model_couplers/coupler_solid_contact.hh
+++ b/src/model/model_couplers/coupler_solid_contact.hh
@@ -1,189 +1,228 @@
/**
* @file coupler_solid_contact_explicit.hh
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Thu Jan 17 2019
* @date last modification: Thu Jan 17 2019
*
* @brief class for coupling of solid mechanics and conatct mechanics
* model in explicit
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "boundary_condition.hh"
#include "contact_mechanics_model.hh"
#include "data_accessor.hh"
#include "model.hh"
#include "solid_mechanics_model.hh"
#include "sparse_matrix.hh"
#include "time_step_solver.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_COUPLER_SOLID_CONTACT_HH__
#define __AKANTU_COUPLER_SOLID_CONTACT_HH__
/* ------------------------------------------------------------------------ */
/* Coupling : Solid Mechanics / Contact Mechanics */
/* ------------------------------------------------------------------------ */
namespace akantu {
template <ElementKind kind, class IntegrationOrderFunctor>
class IntegratorGauss;
template <ElementKind kind> class ShapeLagrange;
class DOFManager;
} // namespace akantu
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
class CouplerSolidContact : public Model,
public DataAccessor<Element>,
public DataAccessor<UInt>,
public BoundaryCondition<CouplerSolidContact> {
/* ------------------------------------------------------------------------ */
/* Constructor/Destructor */
/* ------------------------------------------------------------------------ */
using MyFEEngineType = FEEngineTemplate<IntegratorGauss, ShapeLagrange>;
public:
- CouplerSolidContact(
- SolidMechanicsModel &, ContactMechanicsModel &,
+ CouplerSolidContact(Mesh & mesh,
UInt spatial_dimension = _all_dimensions,
const ID & id = "coupler_solid_contact",
std::shared_ptr<DOFManager> dof_manager = nullptr,
const ModelType model_type = ModelType::_coupler_solid_contact);
~CouplerSolidContact() override;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
protected:
/// initialize completely the model
void initFullImpl(const ModelOptions & options) override;
/// allocate all vectors
void initSolver(TimeStepSolverType, NonLinearSolverType) override;
/// initialize the modelType
void initModel() override;
/// call back for the solver, computes the force residual
void assembleResidual() override;
/// get the type of matrix needed
MatrixType getMatrixType(const ID & matrix_id) override;
/// callback for the solver, this assembles different matrices
void assembleMatrix(const ID & matrix_id) override;
/// callback for the solver, this assembles the stiffness matrix
void assembleLumpedMatrix(const ID & matrix_id) override;
/// get some default values for derived classes
std::tuple<ID, TimeStepSolverType>
getDefaultSolverID(const AnalysisMethod & method) override;
ModelSolverOptions
getDefaultSolverOptions(const TimeStepSolverType & type) const;
/// 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<Element> & elements,
+ const SynchronizationTag & tag) const override;
+
+ void packData(CommunicationBuffer & buffer, const Array<Element> & elements,
+ const SynchronizationTag & tag) const override;
+
+ void unpackData(CommunicationBuffer & buffer, const Array<Element> & elements,
+ const SynchronizationTag & tag) override;
+
+ UInt getNbData(const Array<UInt> & dofs,
+ const SynchronizationTag & tag) const override;
+
+ void packData(CommunicationBuffer & buffer, const Array<UInt> & dofs,
+ const SynchronizationTag & tag) const override;
+
+ void unpackData(CommunicationBuffer & buffer, const Array<UInt> & dofs,
+ const SynchronizationTag & tag) override;
+
+
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
FEEngine & getFEEngineBoundary(const ID & name = "") override;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// return the dimension of the system space
AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt);
/// get the ContactMechanicsModel::displacement vector
AKANTU_GET_MACRO(Displacement, *displacement, Array<Real> &);
/// get the ContactMechanicsModel::increment vector \warn only consistent
/// if ContactMechanicsModel::setIncrementFlagOn has been called before
AKANTU_GET_MACRO(Increment, *displacement_increment, Array<Real> &);
/// get the ContactMechanicsModel::external_force vector (external forces)
AKANTU_GET_MACRO(ExternalForce, *external_force, Array<Real> &);
/// get the ContactMechanicsModel::force vector (external forces)
Array<Real> & getForce() {
AKANTU_DEBUG_WARNING("getForce was maintained for backward compatibility, "
"use getExternalForce instead");
return *external_force;
}
+ /// get the solid mechanics model
+ AKANTU_GET_MACRO(SolidMechanicsModel, *solid, SolidMechanicsModel &);
+
+ /// get the contact mechanics model
+ AKANTU_GET_MACRO(ContactMechanicsModel, *contact, ContactMechanicsModel &);
+
/* ------------------------------------------------------------------------ */
/* 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;
+
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
private:
/// couples external forces between models
void coupleExternalForces();
/// couples stiffness matrices between models
void coupleStiffnessMatrices();
/* ------------------------------------------------------------------------ */
/* Members */
/* ------------------------------------------------------------------------ */
private:
/// solid mechanics model
- SolidMechanicsModel & solid;
+ SolidMechanicsModel * solid{nullptr};
/// contact mechanics model
- ContactMechanicsModel & contact;
+ ContactMechanicsModel * contact{nullptr};
- /// displacements array
+ ///
Array<Real> * displacement{nullptr};
- /// increment of displacement
+ ///
Array<Real> * displacement_increment{nullptr};
-
+
/// external forces array
Array<Real> * external_force{nullptr};
};
} // namespace akantu
#endif /* __COUPLER_SOLID_CONTACT_HH__ */
diff --git a/src/model/model_solver.hh b/src/model/model_solver.hh
index 52724261b..410623aaa 100644
--- a/src/model/model_solver.hh
+++ b/src/model/model_solver.hh
@@ -1,201 +1,201 @@
/**
* @file model_solver.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Wed Jan 31 2018
*
* @brief Class regrouping the common solve interface to the different models
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "integration_scheme.hh"
#include "parsable.hh"
#include "solver_callback.hh"
#include "synchronizer_registry.hh"
/* -------------------------------------------------------------------------- */
#include <set>
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MODEL_SOLVER_HH__
#define __AKANTU_MODEL_SOLVER_HH__
namespace akantu {
class Mesh;
class DOFManager;
class TimeStepSolver;
class NonLinearSolver;
struct ModelSolverOptions;
}
namespace akantu {
class ModelSolver : public Parsable,
public SolverCallback,
public SynchronizerRegistry {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
ModelSolver(Mesh & mesh, const ModelType & type, const ID & id,
- UInt memory_id, std::shared_ptr<DOFManager> dof_manager);
+ UInt memory_id, std::shared_ptr<DOFManager> dof_manager=nullptr);
~ModelSolver() override;
/// initialize the dof manager based on solver type passed in the input file
void initDOFManager();
/// initialize the dof manager based on the used chosen solver type
void initDOFManager(const ID & solver_type);
protected:
/// initialize the dof manager based on the used chosen solver type
void initDOFManager(const ParserSection & section, const ID & solver_type);
/// Callback for the model to instantiate the matricees when needed
virtual void initSolver(TimeStepSolverType /*time_step_solver_type*/,
NonLinearSolverType /*non_linear_solver_type*/) {}
/// get the section in the input file (if it exsits) corresponding to this
/// model
std::tuple<ParserSection, bool> getParserSection();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// solve a step using a given pre instantiated time step solver and
/// non linear solver
virtual void solveStep(const ID & solver_id = "");
/// solve a step using a given pre instantiated time step solver and
/// non linear solver with a user defined callback instead of the
/// model itself /!\ This can mess up everything
virtual void solveStep(SolverCallback & callback,
const ID & solver_id = "");
/// Initialize a time solver that can be used afterwards with its id
void getNewSolver(const ID & solver_id,
TimeStepSolverType time_step_solver_type,
NonLinearSolverType non_linear_solver_type = _nls_auto);
/// set an integration scheme for a given dof and a given solver
void
setIntegrationScheme(const ID & solver_id, const ID & dof_id,
const IntegrationSchemeType & integration_scheme_type,
IntegrationScheme::SolutionType solution_type =
IntegrationScheme::_not_defined);
/// set an externally instantiated integration scheme
void setIntegrationScheme(const ID & solver_id, const ID & dof_id,
IntegrationScheme & integration_scheme,
IntegrationScheme::SolutionType solution_type =
IntegrationScheme::_not_defined);
/* ------------------------------------------------------------------------ */
/* SolverCallback interface */
/* ------------------------------------------------------------------------ */
public:
/// Predictor interface for the callback
void predictor() override;
/// Corrector interface for the callback
void corrector() override;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// Default time step solver to instantiate for this model
virtual TimeStepSolverType getDefaultSolverType() const;
/// Default configurations for a given time step solver
virtual ModelSolverOptions
getDefaultSolverOptions(const TimeStepSolverType & type) const;
/// get access to the internal dof manager
DOFManager & getDOFManager() { return *this->dof_manager; }
/// get the time step of a given solver
Real getTimeStep(const ID & solver_id = "") const;
/// set the time step of a given solver
virtual void setTimeStep(Real time_step, const ID & solver_id = "");
/// set the parameter 'param' of the solver 'solver_id'
// template <typename T>
// void set(const ID & param, const T & value, const ID & solver_id = "");
/// get the parameter 'param' of the solver 'solver_id'
// const Parameter & get(const ID & param, const ID & solver_id = "") const;
/// answer to the question "does the solver exists ?"
bool hasSolver(const ID & solver_id) const;
/// changes the current default solver
void setDefaultSolver(const ID & solver_id);
/// is a default solver defined
bool hasDefaultSolver() const;
/// is an integration scheme set for a given solver and a given dof
bool hasIntegrationScheme(const ID & solver_id, const ID & dof_id) const;
TimeStepSolver & getTimeStepSolver(const ID & solver_id = "");
NonLinearSolver & getNonLinearSolver(const ID & solver_id = "");
const TimeStepSolver & getTimeStepSolver(const ID & solver_id = "") const;
const NonLinearSolver & getNonLinearSolver(const ID & solver_id = "") const;
private:
TimeStepSolver & getSolver(const ID & solver_id);
const TimeStepSolver & getSolver(const ID & solver_id) const;
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
ModelType model_type;
/// Underlying dof_manager (the brain...)
std::shared_ptr<DOFManager> dof_manager;
private:
ID parent_id;
UInt parent_memory_id;
/// Underlying mesh
Mesh & mesh;
/// Default time step solver to use
ID default_solver_id{""};
};
struct ModelSolverOptions {
NonLinearSolverType non_linear_solver_type;
std::map<ID, IntegrationSchemeType> integration_scheme_type;
std::map<ID, IntegrationScheme::SolutionType> solution_type;
};
} // akantu
#endif /* __AKANTU_MODEL_SOLVER_HH__ */
diff --git a/test/test_model/test_contact_mechanics_model/test_coupler/test_contact_coupling.cc b/test/test_model/test_contact_mechanics_model/test_coupler/test_contact_coupling.cc
index 7e859ad0d..f6b5938c9 100644
--- a/test/test_model/test_contact_mechanics_model/test_coupler/test_contact_coupling.cc
+++ b/test/test_model/test_contact_mechanics_model/test_coupler/test_contact_coupling.cc
@@ -1,70 +1,76 @@
/**
* @file test_contact_coupler.cc
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Tue Apr 30 2019
* @date last modification: Tue Apr 30 2019
*
* @brief Test for contact mechanics model 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "coupler_solid_contact.hh"
-#include "non_linear_solver.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
int main(int argc, char *argv[]) {
const UInt spatial_dimension = 2;
initialize("material.dat", argc, argv);
Mesh mesh(spatial_dimension);
mesh.read("coupling.msh");
- SolidMechanicsModel solid(mesh);
+ CouplerSolidContact coupler(mesh);
+
+ auto & solid = coupler.getSolidMechanicsModel();
+ auto & contact = coupler.getContactMechanicsModel();
+
solid.initFull(_analysis_method = _static);
+ contact.initFull(_analysis_method = _implicit_contact);
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "bot_body");
solid.applyBC(BC::Dirichlet::IncrementValue(0.001, _y), "bot_body");
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "top");
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "top");
-
- auto & solver = solid.getNonLinearSolver();
- solver.set("max_iterations", 1);
- solver.set("threshold", 1e-1);
- solver.set("convergence_type", _scc_solution);
- ContactMechanicsModel contact(mesh);
- contact.initFull(_analysis_method = _implicit_contact);
-
- CouplerSolidContact coupler(solid, contact);
coupler.initFull(_analysis_method = _implicit_contact);
+
+ coupler.setBaseName("coupling");
+ coupler.addDumpFieldVector("displacement");
+ coupler.addDumpField("blocked_dofs");
+ coupler.addDumpField("external_force");
+ coupler.addDumpField("internal_force");
+ coupler.addDumpField("grad_u");
+ coupler.addDumpField("stress");
+
coupler.solveStep();
+ contact.dump();
+
return 0;
}