diff --git a/src/model/contact_mechanics/contact_detector.cc b/src/model/contact_mechanics/contact_detector.cc
index 65ff4fd7a..7fc7ef6dd 100644
--- a/src/model/contact_mechanics/contact_detector.cc
+++ b/src/model/contact_mechanics/contact_detector.cc
@@ -1,613 +1,614 @@
/**
* @file contact_detector.cc
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Wed Sep 12 2018
* @date last modification: Fri Sep 21 2018
*
* @brief Mother class for all detection algorithms
*
* @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_detector.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
ContactDetector::ContactDetector(Mesh & mesh, const ID & id,
const MemoryID & memory_id)
: ContactDetector(mesh, mesh.getNodes(), id, memory_id) {}
/* -------------------------------------------------------------------------- */
ContactDetector::ContactDetector(Mesh & mesh, Array<Real> positions,
const ID & id, const MemoryID & memory_id)
: Memory(id, memory_id), Parsable(ParserType::_contact_detector, id),
mesh(mesh), positions(0, mesh.getSpatialDimension()) {
AKANTU_DEBUG_IN();
this->spatial_dimension = mesh.getSpatialDimension();
this->positions.copy(positions);
this->parseSection();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ContactDetector::parseSection() {
const Parser & parser = getStaticParser();
const ParserSection & section =
*(parser.getSubSections(ParserType::_contact_detector).first);
auto type = section.getParameterValue<std::string>("type");
if (type == "implicit") {
this->detection_type = _implicit;
} else if (type == "explicit") {
this->detection_type = _explicit;
} else {
AKANTU_ERROR("Unknown detection type : " << type);
}
}
/* -------------------------------------------------------------------------- */
void ContactDetector::search(Array<ContactElement> & elements,
Array<Real> & gaps, Array<Real> & normals,
Array<Real> & projections) {
UInt surface_dimension = spatial_dimension - 1;
this->mesh.fillNodesToElements(surface_dimension);
this->computeMaximalDetectionDistance();
contact_pairs.clear();
SpatialGrid<UInt> master_grid(spatial_dimension);
SpatialGrid<UInt> slave_grid(spatial_dimension);
this->globalSearch(slave_grid, master_grid);
this->localSearch(slave_grid, master_grid);
this->createContactElements(elements, gaps, normals, projections);
}
/* -------------------------------------------------------------------------- */
void ContactDetector::globalSearch(SpatialGrid<UInt> & slave_grid,
SpatialGrid<UInt> & master_grid) {
auto & master_list = surface_selector->getMasterList();
auto & slave_list = surface_selector->getSlaveList();
BBox bbox_master(spatial_dimension);
this->constructBoundingBox(bbox_master, master_list);
BBox bbox_slave(spatial_dimension);
this->constructBoundingBox(bbox_slave, slave_list);
auto && bbox_intersection = bbox_master.intersection(bbox_slave);
AKANTU_DEBUG_INFO("Intersection BBox " << bbox_intersection);
Vector<Real> center(spatial_dimension);
bbox_intersection.getCenter(center);
Vector<Real> spacing(spatial_dimension);
this->computeCellSpacing(spacing);
master_grid.setCenter(center);
master_grid.setSpacing(spacing);
this->constructGrid(master_grid, bbox_intersection, master_list);
slave_grid.setCenter(center);
slave_grid.setSpacing(spacing);
this->constructGrid(slave_grid, bbox_intersection, slave_list);
// search slave grid nodes in contactelement array and if they exits
// and still have orthogonal projection on its associated master
// facetremove it from the spatial grid or do not consider it for
// local search, maybe better option will be to have spatial grid of
// type node info and one of the variable of node info should be
// facet already exits
// so contact elements will be updated based on the above
// consideration , this means only those contact elements will be
// keep whose slave node is still in intersection bbox and still has
// projection in its master facet
// also if slave node is already exists in contact element and
// orthogonal projection does not exits then search the associated
// master facets with the current master facets within a given
// radius , this is subjected to computational cost as searching
// neighbbor cells can be more effective.
}
/* -------------------------------------------------------------------------- */
void ContactDetector::localSearch(SpatialGrid<UInt> & slave_grid,
SpatialGrid<UInt> & master_grid) {
// local search
// out of these array check each cell for closet node in that cell
// and neighbouring cells find the actual orthogonally closet
// check the projection of slave node on master facets connected to
// the closet master node, if yes update the contact element with
// slave node and master node and master surfaces connected to the
// master node
// these master surfaces will be needed later to update contact
// elements
/// find the closet master node for each slave node
for (auto && cell_id : slave_grid) {
/// loop over all the slave nodes of the current cell
for (auto && slave_node : slave_grid.getCell(cell_id)) {
bool pair_exists = false;
Vector<Real> pos(spatial_dimension);
for (UInt s : arange(spatial_dimension))
pos(s) = this->positions(slave_node, s);
Real closet_distance = std::numeric_limits<Real>::max();
UInt closet_master_node;
/// loop over all the neighboring cells of the current cell
for (auto && neighbor_cell : cell_id.neighbors()) {
/// loop over the data of neighboring cells from master grid
for (auto && master_node : master_grid.getCell(neighbor_cell)) {
/// check for self contact
if (slave_node == master_node)
continue;
bool is_valid = true;
Array<Element> elements;
this->mesh.getAssociatedElements(slave_node, elements);
for (auto & elem : elements) {
if (elem.kind() != _ek_regular)
continue;
Vector<UInt> connectivity =
const_cast<const Mesh &>(this->mesh).getConnectivity(elem);
auto node_iter = std::find(connectivity.begin(), connectivity.end(),
master_node);
if (node_iter != connectivity.end()) {
is_valid = false;
break;
}
}
Vector<Real> pos2(spatial_dimension);
for (UInt s : arange(spatial_dimension))
pos2(s) = this->positions(master_node, s);
Real distance = pos.distance(pos2);
if (distance <= closet_distance and is_valid) {
closet_master_node = master_node;
closet_distance = distance;
pair_exists = true;
}
}
}
if (pair_exists)
contact_pairs.push_back(std::make_pair(slave_node, closet_master_node));
}
}
}
/* -------------------------------------------------------------------------- */
/*void ContactDetector::constructContactMap(
std::map<UInt, ContactElement> & contact_map) {
auto surface_dimension = spatial_dimension - 1;
std::map<UInt, ContactElement> previous_contact_map = contact_map;
contact_map.clear();
auto get_connectivity = [&](auto & slave, auto & master) {
Vector<UInt> master_conn =
const_cast<const Mesh &>(this->mesh).getConnectivity(master);
Vector<UInt> elem_conn(master_conn.size() + 1);
elem_conn[0] = slave;
for (UInt i = 1; i < elem_conn.size(); ++i) {
elem_conn[i] = master_conn[i - 1];
}
return elem_conn;
};
for (auto & pairs : contact_pairs) {
const auto & slave_node = pairs.first;
const auto & master_node = pairs.second;
Array<Element> all_elements;
this->mesh.getAssociatedElements(master_node, all_elements);
Array<Element> boundary_elements;
this->filterBoundaryElements(all_elements, boundary_elements);
Array<Real> gaps(boundary_elements.size(), 1, "gaps");
Array<Real> normals(boundary_elements.size(), spatial_dimension, "normals");
Array<Real> projections(boundary_elements.size(), surface_dimension,
"projections");
auto index = this->computeOrthogonalProjection(
slave_node, boundary_elements, normals, gaps, projections);
if (index == UInt(-1)) {
continue;
}
auto connectivity = get_connectivity(slave_node, boundary_elements[index]);
// assign contact element attributes
contact_map[slave_node].setMaster(boundary_elements[index]);
contact_map[slave_node].setGap(gaps[index]);
contact_map[slave_node].setNormal(
Vector<Real>(normals.begin(spatial_dimension)[index], true));
contact_map[slave_node].setProjection(
Vector<Real>(projections.begin(surface_dimension)[index], true));
contact_map[slave_node].setConnectivity(connectivity);
// tangent computation on master surface
Matrix<Real> tangents(surface_dimension, spatial_dimension);
this->computeTangentsOnElement(contact_map[slave_node].master,
contact_map[slave_node].projection,
contact_map[slave_node].normal, tangents);
contact_map[slave_node].setTangent(tangents);
// friction calculation requires history of previous natural
// projection as well as traction
auto search = previous_contact_map.find(slave_node);
if (search != previous_contact_map.end()) {
auto previous_projection =
previous_contact_map[slave_node].getPreviousProjection();
contact_map[slave_node].setPreviousProjection(previous_projection);
} else {
Vector<Real> previous_projection(surface_dimension, 0.);
contact_map[slave_node].setPreviousProjection(previous_projection);
}
// to ensure the self contact between surface does not lead to
// detection of master element which is closet but not
// orthogonally opposite to the slave surface
bool is_valid_self_contact =
this->checkValidityOfSelfContact(slave_node, contact_map[slave_node]);
if (!is_valid_self_contact) {
contact_map.erase(slave_node);
}
}
contact_pairs.clear();
}*/
/* -------------------------------------------------------------------------- */
void ContactDetector::createContactElements(Array<ContactElement> & contact_elements,
Array<Real> & gaps, Array<Real> & normals,
Array<Real> & projections) {
auto surface_dimension = spatial_dimension - 1;
Real alpha;
switch (detection_type) {
case _explicit: {
alpha = 1.0;
break;
}
case _implicit: {
alpha = -1.0;
break;
}
default:
AKANTU_EXCEPTION(detection_type << " is not a valid contact detection type");
break;
}
for (auto & pairs : contact_pairs) {
const auto & slave_node = pairs.first;
+
Vector<Real> slave(spatial_dimension);
for (UInt s : arange(spatial_dimension))
slave(s) = this->positions(slave_node, s);
const auto & master_node = pairs.second;
Array<Element> elements;
this->mesh.getAssociatedElements(master_node, elements);
auto & gap = gaps.begin()[slave_node];
Vector<Real> normal(normals.begin(spatial_dimension)[slave_node]);
Vector<Real> projection(projections.begin(surface_dimension)[slave_node]);
auto index = GeometryUtils::orthogonalProjection(mesh, positions, slave, elements,
gap, projection, normal, alpha);
// if not a valid projection is found on patch of elements index is -1
if (index == UInt(-1))
continue;
// if not a valid self contact
if (!isValidSelfContact(slave_node, gap, normal))
continue;
// create contact element
contact_elements.push_back(ContactElement(slave_node, elements[index]));
}
contact_pairs.clear();
}
/* -------------------------------------------------------------------------- */
/*UInt ContactDetector::computeOrthogonalProjection(
const UInt & node, const Array<Element> & elements, Array<Real> & normals,
Array<Real> & gaps, Array<Real> & projections) {
Vector<Real> query(spatial_dimension);
for (UInt s : arange(spatial_dimension))
query(s) = this->positions(node, s);
UInt counter = 0;
UInt index = UInt(-1);
Real min_gap = std::numeric_limits<Real>::max();
for (auto && values :
zip(elements, gaps, make_view(normals, spatial_dimension),
make_view(projections, spatial_dimension - 1))) {
const auto & element = std::get<0>(values);
auto & gap = std::get<1>(values);
auto & normal = std::get<2>(values);
auto & projection = std::get<3>(values);
this->computeNormalOnElement(element, normal);
Vector<Real> real_projection(spatial_dimension);
this->computeProjectionOnElement(element, normal, query, projection,
real_projection);
gap = this->computeGap(query, real_projection);
// check if gap is valid or not
// to check this we need normal on master element, vector from
// real projection to slave node, if it is explicit detection
// scheme, we want it to interpenetrate, the dot product should be
// negative and -1.0 and for implciit detection , the dot product
// should be positive and 1.0
bool is_valid = this->checkValidityOfProjection(projection);
auto master_to_slave = query - real_projection;
auto norm = master_to_slave.norm();
if (norm != 0)
master_to_slave /= norm;
Real tolerance = 1e-8;
switch (detection_type) {
case _explicit: {
auto product = master_to_slave.dot(normal);
auto variation = std::abs(product + 1.0);
if (variation <= tolerance and gap <= min_gap and is_valid) {
min_gap = gap;
index = counter;
gap *= -1.0;
}
break;
}
case _implicit: {
auto product = master_to_slave.dot(normal);
auto variation = std::abs(product - 1.0);
if (variation <= tolerance and gap <= min_gap and is_valid) {
min_gap = gap;
index = counter;
gap *= -1.0;
}
}
default:
break;
}
counter++;
}
return index;
}*/
/* -------------------------------------------------------------------------- */
/*void ContactDetector::computeProjectionOnElement(
const Element & element, const Vector<Real> & normal,
const Vector<Real> & query, Vector<Real> & natural_projection,
Vector<Real> & real_projection) {
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type);
Matrix<Real> coords(spatial_dimension, nb_nodes_per_element);
this->coordinatesOfElement(element, coords);
Vector<Real> point(coords(0));
Real alpha = (query - point).dot(normal);
real_projection = query - alpha * normal;
this->computeNaturalProjection(element, real_projection, natural_projection);
}*/
/* -------------------------------------------------------------------------- */
/*void ContactDetector::computeNaturalProjection(
const Element & element, Vector<Real> & real_projection,
Vector<Real> & natural_projection) {
const ElementType & type = element.type;
UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type);
UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage();
Matrix<Real> nodes_coord(spatial_dimension, nb_nodes_per_element);
mesh.extractNodalValuesFromElement(this->positions, nodes_coord.storage(),
elem_val +
element.element * nb_nodes_per_element,
nb_nodes_per_element, spatial_dimension);
#define GET_NATURAL_COORDINATE(type) \
ElementClass<type>::inverseMap(real_projection, nodes_coord, \
natural_projection)
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NATURAL_COORDINATE);
#undef GET_NATURAL_COORDINATE
}*/
/* -------------------------------------------------------------------------- */
/*void ContactDetector::computeTangentsOnElement(const Element & el,
Vector<Real> & projection,
Matrix<Real> & tangents) {
const ElementType & type = el.type;
UInt nb_nodes_master = Mesh::getNbNodesPerElement(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> coords(spatial_dimension, nb_nodes_master);
coordinatesOfElement(el, coords);
tangents.mul<false, true>(shapes_derivatives, coords);
auto temp_tangents = tangents.transpose();
for (UInt i = 0; i < spatial_dimension - 1; ++i) {
auto temp = Vector<Real>(temp_tangents(i));
temp_tangents(i) = temp.normalize();
}
tangents = temp_tangents.transpose();
}*/
/* -------------------------------------------------------------------------- */
/*void ContactDetector::computeTangentsOnElement(const Element & el,
Vector<Real> & projection,
Vector<Real> & normal,
Matrix<Real> & tangents) {
const ElementType & type = el.type;
UInt nb_nodes_master = Mesh::getNbNodesPerElement(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> coords(spatial_dimension, nb_nodes_master);
coordinatesOfElement(el, coords);
tangents.mul<false, true>(shapes_derivatives, coords);
auto temp_tangents = tangents.transpose();
for (UInt i = 0; i < spatial_dimension - 1; ++i) {
auto temp = Vector<Real>(temp_tangents(i));
temp_tangents(i) = temp.normalize();
}
tangents = temp_tangents.transpose();
// to ensure that direction of tangents are correct, cross product
// of tangents should give the normal vector computed earlier
switch (spatial_dimension) {
case 2: {
Vector<Real> e_z(3);
e_z[0] = 0.;
e_z[1] = 0.;
e_z[2] = 1.;
Vector<Real> tangent(3);
tangent[0] = tangents(0, 0);
tangent[1] = tangents(0, 1);
tangent[2] = 0.;
auto exp_normal = e_z.crossProduct(tangent);
auto & cal_normal = normal;
auto ddot = cal_normal.dot(exp_normal);
if (ddot < 0) {
tangents *= -1.0;
}
break;
}
case 3: {
auto tang_trans = tangents.transpose();
auto tang1 = Vector<Real>(tang_trans(0));
auto tang2 = Vector<Real>(tang_trans(1));
auto tang1_cross_tang2 = tang1.crossProduct(tang2);
auto exp_normal = tang1_cross_tang2 / tang1_cross_tang2.norm();
auto & cal_normal = normal;
auto ddot = cal_normal.dot(exp_normal);
if (ddot < 0) {
tang_trans(1) *= -1.0;
}
tangents = tang_trans.transpose();
break;
}
default:
break;
}
}*/
} // namespace akantu
diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc
index 60eb838d6..522c5d18a 100644
--- a/src/model/contact_mechanics/contact_mechanics_model.cc
+++ b/src/model/contact_mechanics/contact_mechanics_model.cc
@@ -1,682 +1,685 @@
/**
* @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"
#include "group_manager_inline_impl.hh"
#ifdef AKANTU_USE_IOHELPER
#include "dumper_iohelper_paraview.hh"
#endif
/* -------------------------------------------------------------------------- */
#include <algorithm>
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
ContactMechanicsModel::ContactMechanicsModel(
Mesh & mesh, 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) {
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, id + ":contact_detector");
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
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) {
// for alloc type of solvers
this->allocNodalField(this->displacement, spatial_dimension, "displacement");
this->allocNodalField(this->displacement_increment, spatial_dimension,
"displacement_increment");
this->allocNodalField(this->internal_force, spatial_dimension,
"internal_force");
this->allocNodalField(this->external_force, spatial_dimension,
"external_force");
this->allocNodalField(this->normal_force, spatial_dimension,
"normal_force");
this->allocNodalField(this->tangential_force, spatial_dimension,
"tangential_force");
this->allocNodalField(this->gaps, 1, "gaps");
this->allocNodalField(this->nodal_area, 1, "areas");
this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs");
this->allocNodalField(this->stick_or_slip, 1, "stick_or_slip");
this->allocNodalField(this->previous_master_elements, 1, "previous_master_elements");
this->allocNodalField(this->normals, spatial_dimension, "normals");
auto surface_dimension = spatial_dimension - 1;
this->allocNodalField(this->tangents, surface_dimension*spatial_dimension,
"tangents");
this->allocNodalField(this->previous_tangents, surface_dimension*spatial_dimension,
"previous_tangents");
this->allocNodalField(this->projections, surface_dimension,
"projections");
this->allocNodalField(this->previous_projections, surface_dimension,
"previous_projections");
this->allocNodalField(this->stick_projections, surface_dimension,
"stick_projections");
this->allocNodalField(this->tangential_tractions, surface_dimension,
"tangential_tractions");
this->allocNodalField(this->previous_tangential_tractions, surface_dimension,
"previous_tangential_tractions");
// todo register multipliers as dofs for lagrange multipliers
}
/* -------------------------------------------------------------------------- */
std::tuple<ID, TimeStepSolverType>
ContactMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) {
switch (method) {
case _explicit_lumped_mass: {
return std::make_tuple("explicit_lumped",
TimeStepSolverType::_dynamic_lumped);
}
case _explicit_consistent_mass: {
return std::make_tuple("explicit", TimeStepSolverType::_dynamic);
}
case _static: {
return std::make_tuple("static", TimeStepSolverType::_static);
}
case _implicit_dynamic: {
return std::make_tuple("implicit", TimeStepSolverType::_dynamic);
}
default:
return std::make_tuple("unknown", TimeStepSolverType::_not_defined);
}
}
/* -------------------------------------------------------------------------- */
ModelSolverOptions ContactMechanicsModel::getDefaultSolverOptions(
const TimeStepSolverType & type) const {
ModelSolverOptions options;
switch (type) {
case TimeStepSolverType::_dynamic: {
options.non_linear_solver_type = NonLinearSolverType::_lumped;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_central_difference;
options.solution_type["displacement"] = IntegrationScheme::_acceleration;
break;
}
case TimeStepSolverType::_dynamic_lumped: {
options.non_linear_solver_type = NonLinearSolverType::_lumped;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_central_difference;
options.solution_type["displacement"] = IntegrationScheme::_acceleration;
break;
}
case TimeStepSolverType::_static: {
options.non_linear_solver_type =
NonLinearSolverType::_newton_raphson_contact;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_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");
UInt nb_nodes = mesh.getNbNodes();
this->internal_force->clear();
this->normal_force->clear();
this->tangential_force->clear();
internal_force->resize(nb_nodes, 0.);
normal_force->resize(nb_nodes, 0.);
tangential_force->resize(nb_nodes, 0.);
// assemble the forces due to contact
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);
+ //assemble(_ghost);
+ // TODO : uncomment when developing code for parallelization,
+ // currently it addes the force twice for not ghost elements
+ // hence source of error
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::search() {
// save the previous state
this->savePreviousState();
contact_elements.clear();
contact_elements.resize(0);
// this needs to be resized if cohesive elements are added
UInt nb_nodes = mesh.getNbNodes();
auto resize_arrays = [&](auto & internal_array) {
internal_array->clear();
internal_array->resize(nb_nodes, 0.);
};
resize_arrays(gaps);
resize_arrays(normals);
resize_arrays(tangents);
resize_arrays(projections);
resize_arrays(tangential_tractions);
resize_arrays(stick_or_slip);
this->detector->search(contact_elements, *gaps,
*normals, *projections);
for (auto & gap : *gaps) {
if (gap < 0)
gap = std::abs(gap);
else
gap = -gap;
}
if (contact_elements.size() != 0) {
this->computeNodalAreas();
}
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::savePreviousState() {
AKANTU_DEBUG_IN();
// saving previous natural projections
previous_projections->clear();
previous_projections->resize(projections->size(), 0.);
(*previous_projections).copy(*projections);
// saving previous tangents
previous_tangents->clear();
previous_tangents->resize(tangents->size(), 0.);
(*previous_tangents).copy(*tangents);
// saving previous tangential traction
previous_tangential_tractions->clear();
previous_tangential_tractions->resize(tangential_tractions->size(), 0.);
(*previous_tangential_tractions).copy(*tangential_tractions);
previous_master_elements->clear();
previous_master_elements->resize(projections->size());
for (auto & element : contact_elements) {
(*previous_master_elements)[element.slave] = element.master;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::computeNodalAreas() {
UInt nb_nodes = mesh.getNbNodes();
nodal_area->clear();
external_force->clear();
nodal_area->resize(nb_nodes, 0.);
external_force->resize(nb_nodes, 0.);
auto & fem_boundary = this->getFEEngineBoundary();
fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost);
fem_boundary.computeNormalsOnIntegrationPoints(_ghost);
// TODO find a better method to compute nodal area
switch (spatial_dimension) {
case 1: {
for (auto && area : *nodal_area) {
area = 1.;
}
break;
}
case 2:
case 3: {
this->applyBC(
BC::Neumann::FromHigherDim(Matrix<Real>::eye(spatial_dimension, 1)),
mesh.getElementGroup("contact_surface"));
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);
}
break;
}
default:
break;
}
this->external_force->clear();
}
/* -------------------------------------------------------------------------- */
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"));
}
for (auto & resolution : resolutions) {
resolution->assembleStiffnessMatrix(_not_ghost);
}
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::assembleLumpedMatrix(const ID & /*matrix_id*/) {
AKANTU_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::beforeSolveStep() {
for (auto & resolution : resolutions)
resolution->beforeSolveStep();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::afterSolveStep(bool converged) {
for (auto & resolution : resolutions)
resolution->afterSolveStep(converged);
}
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
/* -------------------------------------------------------------------------- */
std::shared_ptr<dumpers::Field>
ContactMechanicsModel::createNodalFieldBool(const std::string &,
const std::string &, bool) {
return nullptr;
}
/* -------------------------------------------------------------------------- */
std::shared_ptr<dumpers::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->internal_force;
real_nodal_fields["normal_force"] = this->normal_force;
real_nodal_fields["tangential_force"] = this->tangential_force;
real_nodal_fields["blocked_dofs"] = this->blocked_dofs;
real_nodal_fields["normals"] = this->normals;
real_nodal_fields["tangents"] = this->tangents;
real_nodal_fields["gaps"] = this->gaps;
real_nodal_fields["areas"] = this->nodal_area;
real_nodal_fields["stick_or_slip"] = this->stick_or_slip;
std::shared_ptr<dumpers::Field> field;
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
/* -------------------------------------------------------------------------- */
std::shared_ptr<dumpers::Field>
ContactMechanicsModel::createNodalFieldBool(const std::string &,
const std::string &, bool) {
return nullptr;
}
/* -------------------------------------------------------------------------- */
std::shared_ptr<dumpers::Field>
ContactMechanicsModel::createNodalFieldReal(const std::string & ,
const std::string & , bool) {
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;
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/resolutions/resolution_penalty.cc b/src/model/contact_mechanics/resolutions/resolution_penalty.cc
index 15e58e280..4594fa08a 100644
--- a/src/model/contact_mechanics/resolutions/resolution_penalty.cc
+++ b/src/model/contact_mechanics/resolutions/resolution_penalty.cc
@@ -1,695 +1,695 @@
/**
* @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_n", epsilon_n, 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(const ContactElement & element,
Vector<Real> & force) {
force.clear();
auto & gaps = model.getGaps();
auto & projections = model.getProjections();
auto & normals = model.getNormals();
auto surface_dimension = spatial_dimension - 1;
Real gap(gaps.begin()[element.slave]);
Vector<Real> normal(normals.begin(spatial_dimension)[element.slave]);
Vector<Real> projection(projections.begin(surface_dimension)[element.slave]);
// compute normal traction
Real p_n = macaulay(gap) * epsilon_n;
// compute first variation of gap
auto nb_nodes = element.getNbNodes();
Vector<Real> delta_gap(nb_nodes * spatial_dimension);
ResolutionUtils::firstVariationNormalGap(element, projection, normal, delta_gap);
// compute normal force
auto & nodal_area = const_cast<Array<Real> &>(model.getNodalArea());
for (UInt i : arange(force.size()))
force[i] += delta_gap[i] * p_n * nodal_area[element.slave];
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeTangentialForce(const ContactElement & element,
Vector<Real> & force) {
if (mu == 0)
return;
UInt surface_dimension = spatial_dimension - 1;
// compute tangents
auto & projections = model.getProjections();
Vector<Real> projection(projections.begin(surface_dimension)[element.slave]);
auto & normals = model.getNormals();
Vector<Real> normal(normals.begin(spatial_dimension)[element.slave]);
auto & tangents = model.getTangents();
Matrix<Real> covariant_basis(tangents.begin(surface_dimension,
spatial_dimension)[element.slave]);
GeometryUtils::covariantBasis(model.getMesh(), model.getContactDetector().getPositions(),
element.master, normal, projection, covariant_basis);
// check for no-contact to contact condition
auto & previous_master_elements = model.getPreviousMasterElements();
auto & previous_element = previous_master_elements[element.slave];
if (previous_element.type == _not_defined)
return;
// compute tangential traction using return map algorithm
auto & tangential_tractions = model.getTangentialTractions();
Vector<Real> tangential_traction(tangential_tractions.begin(surface_dimension)[element.slave]);
this->computeTangentialTraction(element, covariant_basis,
tangential_traction);
// compute first variation of natural coordinate
auto & gaps = model.getGaps();
auto & gap = gaps.begin()[element.slave];
auto nb_nodes = element.getNbNodes();
Array<Real> delta_xi(nb_nodes * spatial_dimension, surface_dimension);
ResolutionUtils::firstVariationNaturalCoordinate(element, covariant_basis,
projection, normal, gap, delta_xi);
// compute tangential force
auto & nodal_area = const_cast<Array<Real> &>(model.getNodalArea());
for (auto && values : zip(tangential_traction,
make_view(delta_xi, delta_xi.size()))) {
auto & traction_alpha = std::get<0>(values);
auto & delta_xi_alpha = std::get<1>(values);
force += delta_xi_alpha * traction_alpha * nodal_area[element.slave];
}
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeTangentialTraction(const ContactElement & element,
const Matrix<Real> & covariant_basis,
Vector<Real> & traction_tangential) {
UInt surface_dimension = spatial_dimension - 1;
auto & gaps = model.getGaps();
auto & gap = gaps.begin()[element.slave];
// Return map algorithm is employed
// compute trial traction
Vector<Real> traction_trial(surface_dimension);
this->computeTrialTangentialTraction(element, covariant_basis, traction_trial);
// compute norm of trial traction
Real traction_trial_norm = 0;
auto inv_A = GeometryUtils::contravariantMetricTensor(covariant_basis);
for (auto i : arange(surface_dimension)) {
for (auto j : arange(surface_dimension)) {
traction_trial_norm += traction_trial[i] * traction_trial[j] * inv_A(i, j);
}
}
traction_trial_norm = sqrt(traction_trial_norm);
// check stick or slip condition
auto & stick_or_slip = model.getStickSlip();
auto & cond = stick_or_slip.begin()[element.slave];
Real p_n = macaulay(gap) * epsilon_n;
bool stick = (traction_trial_norm <= mu * p_n) ? true : false;
if (stick) {
cond = 1;
computeStickTangentialTraction(element, traction_trial, traction_tangential);
} else {
cond = 2;
computeSlipTangentialTraction(element, covariant_basis, traction_trial,
traction_tangential);
}
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeTrialTangentialTraction(const ContactElement & element,
const Matrix<Real> & covariant_basis,
Vector<Real> & traction) {
UInt surface_dimension = spatial_dimension - 1;
/*auto & projections = model.getProjections();
Vector<Real> contravariant_projection(projections.begin(surface_dimension)[element.slave]);
auto & stick_projections = model.getStickProjections();
Vector<Real> covariant_stick(stick_projections.begin(surface_dimension)[element.slave]);
Matrix<Real> contravariant_basis(surface_dimension, spatial_dimension);
GeometryUtils::contravariantBasis(covariant_basis, contravariant_basis);
Vector<Real> covariant_projection(surface_dimension);
for (auto && values : zip(covariant_projection,
contravariant_projection,
contravariant_basis.transpose())) {
auto & temp = std::get<0>(values);
Vector<Real> contravariant_alpha(std::get<2>(values));
temp = contravariant_alpha.dot(contravariant_alpha);
temp *= std::get<1>(values);
}
auto covariant_slip = covariant_projection - covariant_stick;
traction.mul<true>(contravariant_basis, covariant_slip, epsilon_t);*/
auto & projections = model.getProjections();
Vector<Real> current_projection(projections.begin(surface_dimension)[element.slave]);
auto & previous_projections = model.getPreviousProjections();
Vector<Real> previous_projection(previous_projections.begin(surface_dimension)[element.slave]);
// method from Laursen et. al.
/*auto covariant_metric_tensor = GeometryUtils::covariantMetricTensor(covariant_basis);
auto increment_projection = current_projection - previous_projection;
traction.mul<false>(covariant_metric_tensor, increment_projection, epsilon_t);
auto & previous_tangential_tractions = model.getPreviousTangentialTractions();
Vector<Real> previous_traction(previous_tangential_tractions.begin(surface_dimension)[element.slave]);
traction = previous_traction + traction;*/
// method from Schweizerhof
auto covariant_metric_tensor = GeometryUtils::covariantMetricTensor(covariant_basis);
auto & previous_tangential_tractions = model.getPreviousTangentialTractions();
Vector<Real> previous_traction(previous_tangential_tractions.begin(surface_dimension)[element.slave]);
auto & previous_tangents = model.getPreviousTangents();
Matrix<Real> previous_covariant_basis(previous_tangents.begin(surface_dimension,
spatial_dimension)[element.slave]);
auto previous_contravariant_metric_tensor =
GeometryUtils::contravariantMetricTensor(previous_covariant_basis);
auto current_tangent = covariant_basis.transpose();
auto previous_tangent = previous_covariant_basis.transpose();
for (auto alpha :arange(surface_dimension)) {
Vector<Real> tangent_alpha(current_tangent(alpha));
for (auto gamma : arange(surface_dimension)) {
for (auto beta : arange(surface_dimension)) {
Vector<Real> tangent_beta(previous_tangent(beta));
auto t_alpha_t_beta = tangent_beta.dot(tangent_alpha);
traction[alpha] += previous_traction[gamma]*previous_contravariant_metric_tensor(gamma, beta)*t_alpha_t_beta;
}
}
}
auto & previous_master_elements = model.getPreviousMasterElements();
auto & previous_element = previous_master_elements[element.slave];
Vector<Real> previous_real_projection(spatial_dimension);
GeometryUtils::realProjection(model.getMesh(), model.getContactDetector().getPositions(),
previous_element, previous_projection,
previous_real_projection);
Vector<Real> current_real_projection(spatial_dimension);
GeometryUtils::realProjection(model.getMesh(), model.getContactDetector().getPositions(),
element.master, current_projection,
current_real_projection);
auto increment_real = current_real_projection - previous_real_projection;
Vector<Real> increment_xi(surface_dimension);
auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis);
// increment in natural coordinate
for (auto beta : arange(surface_dimension)) {
for (auto gamma: arange(surface_dimension)) {
auto temp = increment_real.dot(current_tangent(gamma));
temp *= contravariant_metric_tensor(beta, gamma);
increment_xi[beta] += temp;
}
}
Vector<Real> temp(surface_dimension);
temp.mul<false>(covariant_metric_tensor, increment_xi, epsilon_t);
traction -= temp;
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeStickTangentialTraction(const ContactElement & /*element*/,
Vector<Real> & traction_trial,
Vector<Real> & traction_tangential) {
traction_tangential = traction_trial;
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeSlipTangentialTraction(const ContactElement & element,
const Matrix<Real> & covariant_basis,
Vector<Real> & traction_trial,
Vector<Real> & traction_tangential) {
UInt surface_dimension = spatial_dimension - 1;
auto & gaps = model.getGaps();
auto & gap = gaps.begin()[element.slave];
// compute norm of trial traction
Real traction_trial_norm = 0;
auto inv_A = GeometryUtils::contravariantMetricTensor(covariant_basis);
for (auto i : arange(surface_dimension)) {
for (auto j : arange(surface_dimension)) {
traction_trial_norm += traction_trial[i] * traction_trial[j] * inv_A(i, j);
}
}
traction_trial_norm = sqrt(traction_trial_norm);
auto slip_direction = traction_trial;
slip_direction /= traction_trial_norm;
Real p_n = epsilon_n * macaulay(gap);
traction_tangential = slip_direction;
traction_tangential *= mu*p_n;
/*
auto & stick_projections = model.getStickProjections();
Vector<Real> covariant_stick(stick_projections.begin(surface_dimension)[element.slave]);
auto slip = macaulay(traction_trial_norm - mu * p_n);
slip /= epsilon_t;
Vector<Real> covariant_slip(surface_dimension);
for (auto && values : zip(covariant_slip,
covariant_basis.transpose())) {
auto & slip_alpha = std::get<0>(values);
Vector<Real> covariant_alpha(std::get<1>(values));
slip_alpha = slip_direction.dot(covariant_alpha);
}
covariant_stick += slip * covariant_slip;*/
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeNormalModuli(const ContactElement & element,
Matrix<Real> & stiffness) {
auto surface_dimension = spatial_dimension - 1;
auto & gaps = model.getGaps();
Real gap(gaps.begin()[element.slave]);
auto & projections = model.getProjections();
Vector<Real> projection(projections.begin(surface_dimension)[element.slave]);
auto & nodal_areas = model.getNodalArea();
auto & nodal_area = nodal_areas.begin()[element.slave];
auto & normals = model.getNormals();
Vector<Real> normal(normals.begin(spatial_dimension)[element.slave]);
auto & tangents = model.getTangents();
Matrix<Real> covariant_basis(tangents.begin(surface_dimension,
spatial_dimension)[element.slave]);
auto & mesh = model.getMesh();
// construct A matrix
const ElementType & type = element.master.type;
UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type);
Vector<Real> shapes(nb_nodes_per_element);
#define GET_SHAPE_NATURAL(type) \
ElementClass<type>::computeShapes(projection, shapes)
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL);
#undef GET_SHAPE_NATURAL
UInt nb_nodes_per_contact = element.getNbNodes();
Matrix<Real> A(spatial_dimension, spatial_dimension*nb_nodes_per_contact);
for (auto i : arange(nb_nodes_per_contact)) {
for (auto j : arange(spatial_dimension)) {
if (i == 0) {
A(j, i*spatial_dimension + j) = 1;
continue;
}
A(j, i*spatial_dimension + j) = -shapes[i-1];
}
}
// contruct the main part of normal matrix
Matrix<Real> k_main(nb_nodes_per_contact*spatial_dimension, nb_nodes_per_contact*spatial_dimension);
Matrix<Real> n_outer_n(spatial_dimension, spatial_dimension);
Matrix<Real> mat_n(normal.storage(), normal.size(), 1.);
n_outer_n.mul<false, true>(mat_n, mat_n);
Matrix<Real> tmp(spatial_dimension, spatial_dimension*nb_nodes_per_contact);
tmp.mul<false, false>(n_outer_n, A);
k_main.mul<true, false>(A, tmp);
- k_main *= epsilon_n * heaviside(gap) * nodal_area;
+ k_main *= epsilon_n * heaviside(gap); // * nodal_area;
stiffness += k_main;
}
/* -------------------------------------------------------------------------- */
/*void ResolutionPenalty::computeNormalModuli(const ContactElement & element,
Matrix<Real> & ddelta_g,
Vector<Real> & delta_g, Matrix<Real> & stiffness) {
// from Vlad
auto & gaps = model.getGaps();
auto & gap = gaps.begin()[element.slave];
auto & nodal_areas = model.getNodalArea();
auto & nodal_area = nodal_areas.begin()[element.slave];
Matrix<Real> tmp(delta_g.storage(), delta_g.size(), 1);
Matrix<Real> mat_delta_g(delta_g.size(), delta_g.size());
Real heaviside_gap = heaviside(gap);
mat_delta_g.mul<false, true>(tmp, tmp, heaviside_gap);
Real macaulay_gap = macaulay(gap);
ddelta_g *= macaulay_gap;
stiffness += mat_delta_g + ddelta_g;
stiffness *= epsilon_n * nodal_area;
}*/
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeTangentialModuli(const ContactElement & /*element*/,
Matrix<Real> & /*ddelta_g*/,
Vector<Real> & /*delta_g*/,
Matrix<Real> & /*stiffness*/){
if (mu == 0)
return;
}
/* -------------------------------------------------------------------------- */
/*void ResolutionPenalty::computeNormalModuli(Matrix<Real> & ke,
Array<Real> & n_alpha,
Array<Real> & d_alpha,
Vector<Real> & n,
ContactElement & element) {
auto & gaps = model.getGaps();
auto & gap = gaps.begin()[element.slave];
Real tn = gap * epsilon_n;
tn = macaulay(tn);
Matrix<Real> n_mat(n.storage(), n.size(), 1);
ke.mul<false, true>(n_mat, n_mat);
ke *= epsilon_n * 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());
tmp2.mul<false, true>(ds_mat, ns_mat);
ke -= (tmp1 + tmp2) * tn;
}
auto surface_dimension = spatial_dimension - 1;
Matrix<Real> m_alpha_beta(surface_dimension, surface_dimension);
ResolutionUtils::computeMetricTensor(m_alpha_beta, element.tangents);
for (auto && values :
zip(arange(surface_dimension), make_view(n_alpha, n_alpha.size()))) {
auto & s = std::get<0>(values);
auto & n_s = std::get<1>(values);
Matrix<Real> ns_mat(n_s.storage(), n_s.size(), 1);
for (auto && tuple :
zip(arange(surface_dimension), make_view(n_alpha, n_alpha.size()))) {
auto & t = std::get<0>(tuple);
auto & n_t = std::get<1>(tuple);
Matrix<Real> nt_mat(n_t.storage(), n_t.size(), 1);
Matrix<Real> tmp3(n_s.size(), n_s.size());
tmp3.mul<false, true>(ns_mat, nt_mat);
tmp3 *= m_alpha_beta(s, t) * tn * element.gap;
ke += tmp3;
}
}
}*/
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeFrictionalModuli(
Matrix<Real> & /*ke*/, Array<Real> & t_alpha_beta,
Array<Real> & n_alpha_beta, Array<Real> & /*n_alpha*/,
Array<Real> & d_alpha, Matrix<Real> & phi, Vector<Real> & n,
ContactElement & element) {
/*auto k_common =
computeCommonModuli(t_alpha_beta, n_alpha_beta, d_alpha, n, element);
const auto & type = element.master.type;
const auto & conn = element.connectivity;
auto surface_dimension = Mesh::getSpatialDimension(type);
auto spatial_dimension = surface_dimension + 1;
Matrix<Real> m_alpha_beta(surface_dimension, surface_dimension);
ResolutionUtils::computeMetricTensor(m_alpha_beta, element.tangents);
Array<Real> g_alpha(conn.size() * spatial_dimension, surface_dimension);
ResolutionUtils::computeGalpha(g_alpha, t_alpha_beta, d_alpha, phi, element);*/
/*Matrix<Real> k_t;
bool stick = computeFrictionalTraction(m_alpha_beta, element);
if(stick)
k_t = computeStickModuli(g_alpha, d_alpha, m_alpha_beta);
else
k_t = computeSlipModuli(g_alpha, d_alpha, m_alpha_beta, element);*/
}
/* -------------------------------------------------------------------------- */
Array<Real> ResolutionPenalty::computeCommonModuli(
Array<Real> & /*t_alpha_beta*/, Array<Real> & /*n_alpha_beta*/,
Array<Real> & d_alpha, Vector<Real> & /*n*/, ContactElement & /*element*/) {
Array<Real> kt_alpha(spatial_dimension - 1, d_alpha.size() * d_alpha.size(),
"k_T_alpha");
// auto t_alpha_beta_size = t_alpha_beta.size() * (spatial_dimension - 1);
// auto & tangents = element.tangents;
/*for(auto && values :
zip(tangents.transpose(),
make_view(kt_alpha, kt_alpha.size()),
make_view(t_alpha_beta, t_alpha_beta_size),
make_view(n_alpha_beta, n_alpha_beta_size))) {
auto & tangent_s = std::get<0>(values);
auto & kt_s = std::get<1>(values);
auto & t_alpha_s = std::get<2>(values);
auto & n_alpha_s = std::get<3>(values);
Matrix<Real> kt_s_mat(kt_s.storage(), d_alpha.size(), d_alpha.size());
// loop over beta
for(auto && tuple :
zip(make_view(d_alpha, d_alpha.size()),
make_view(n_alpha_ ))) {
auto & d_s = std::get<0>(tuple);
Matrix<Real> tmp(d_s.size(), d_s.size());
// loop over gamma
for(auto && entry :
make_view(d_alpha, d_alpha.size())) {
auto & d_t = std::get<0>(entry);
// compute constant
Matrix<Real> tmp2(d_t.size(), d_t.size());
tmp2.mul<false, true>(d_s, d_t);
kt_s_mat += tmp2;
}
}
}*/
return kt_alpha;
}
/* -------------------------------------------------------------------------- */
Matrix<Real> ResolutionPenalty::computeStickModuli(
Array<Real> & g_alpha, Array<Real> & d_alpha, Matrix<Real> & m_alpha_beta) {
/*Matrix<Real> k_stick(d_alpha.size(), d_alpha.size());
for (auto && values :
zip(arange(d_alpha.getNbComponent()), make_view(d_alpha, d_alpha.size()),
make_view(g_alpha, g_alpha.size()))) {
auto & s = std::get<0>(values);
auto & d_s = std::get<1>(values);
auto & g_s = std::get<2>(values);
Matrix<Real> ds_mat(d_s.storage(), d_s.size(), 1);
Matrix<Real> gs_mat(g_s.storage(), g_s.size(), 1);
Matrix<Real> tmp1(d_s.size(), d_s.size());
tmp1.mul<false, true>(ds_mat, gs_mat);
k_stick += tmp1;
for (auto && tuple : enumerate(make_view(d_alpha, d_alpha.size()))) {
auto & t = std::get<0>(tuple);
auto & d_t = std::get<1>(tuple);
Matrix<Real> dt_mat(d_t.storage(), d_t.size(), 1);
Matrix<Real> tmp2(d_t.size(), d_t.size());
tmp2.mul<false, true>(ds_mat, dt_mat);
k_stick += tmp2 * m_alpha_beta(s, t);
}
}
k_stick *= epsilon_t;
return k_stick;*/
}
/* -------------------------------------------------------------------------- */
Matrix<Real> ResolutionPenalty::computeSlipModuli(
Array<Real> & /*g_alpha*/, Array<Real> & d_alpha,
Matrix<Real> & /*m_alpha_beta*/, ContactElement & /*element*/) {
/*Real tn = element.gap * epsilon;
tn = macaulay(tn);
Real factor;
factor = epsilon_t * mu * tn;
auto p_t = element.traction;
p_t /= p_t.norm();*/
Matrix<Real> k_slip(d_alpha.size(), d_alpha.size());
/*
// loop for alpha
for(auto && value :
make_view(d_alpha, d_alpha.size())) {
auto & d_s = std::get<0>(value);
// loop for beta
for(auto && tuple :
zip(arange(spatial_dimension - 1),
make_view(d_alpha, d_alpha.size()),
make_view(g_alpha, g_alpha.size()))) {
auto & beta = std::get<0>(tuple);
auto & d_beta = std::get<1>(tuple);
auto & g_beta = std::get<2>(tuple);
// loop for gamma
for(auto && entry :
zip(arange(spatial_dimension - 1),
make_view(d_alpha, d_alpha.size()))) {
auto & gamma = std::get<0>(entry);
auto & d_gamma = std::get<1>(entry);
}
}
}*/
return k_slip;
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::beforeSolveStep() {
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::afterSolveStep(bool converged) {
- auto method = model.getAnalysisMethod();
+ /*auto method = model.getAnalysisMethod();
if (method == _explicit_lumped_mass) {
return ;
}
auto & K =
const_cast<SparseMatrix &>(model.getDOFManager().getMatrix("K"));
auto k_min = K.min();
auto roundoff_error = 1e-17;
const auto blocked_dofs = model.getDOFManager().getBlockedDOFs("displacement");
Real nb_unknowns = 0;
for (auto & bld : make_view(blocked_dofs)) {
if (not bld)
nb_unknowns++;
}
auto max_epsilon_n = k_min / sqrt(nb_unknowns * roundoff_error);
if (epsilon_n > max_epsilon_n)
- epsilon_n = max_epsilon_n;
+ epsilon_n = max_epsilon_n;*/
}
INSTANTIATE_RESOLUTION(penalty, ResolutionPenalty);
} // namespace akantu
diff --git a/src/model/model_couplers/coupler_solid_contact.cc b/src/model/model_couplers/coupler_solid_contact.cc
index a2fe4a987..120495f2e 100644
--- a/src/model/model_couplers/coupler_solid_contact.cc
+++ b/src/model/model_couplers/coupler_solid_contact.cc
@@ -1,535 +1,534 @@
/**
* @file coupler_solid_contact.cc
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Thu Jan 17 2019
* @date last modification: Thu May 22 2019
*
* @brief class for coupling of solid mechanics and conatct 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 "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(
Mesh & mesh, 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) {
AKANTU_DEBUG_IN();
this->registerFEEngineObject<MyFEEngineType>("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", 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);
solid->initFull( _analysis_method = this->method);
contact->initFull(_analysis_method = this->method);
}
/* -------------------------------------------------------------------------- */
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 time_step_solver_type,
NonLinearSolverType non_linear_solver_type) {
auto & solid_model_solver =
aka::as_type<ModelSolver>(*solid);
solid_model_solver.initSolver(time_step_solver_type, non_linear_solver_type);
auto & contact_model_solver =
aka::as_type<ModelSolver>(*contact);
contact_model_solver.initSolver(time_step_solver_type, non_linear_solver_type);
}
/* -------------------------------------------------------------------------- */
std::tuple<ID, TimeStepSolverType>
CouplerSolidContact::getDefaultSolverID(const AnalysisMethod & method) {
switch (method) {
case _explicit_lumped_mass: {
return std::make_tuple("explicit_lumped",
TimeStepSolverType::_dynamic_lumped);
}
case _explicit_consistent_mass: {
return std::make_tuple("explicit", TimeStepSolverType::_dynamic);
}
case _static: {
return std::make_tuple("static", TimeStepSolverType::_static);
}
case _implicit_dynamic: {
return std::make_tuple("implicit", TimeStepSolverType::_dynamic);
}
default:
return std::make_tuple("unknown", TimeStepSolverType::_not_defined);
}
}
/* -------------------------------------------------------------------------- */
TimeStepSolverType CouplerSolidContact::getDefaultSolverType() const {
return TimeStepSolverType::_dynamic_lumped;
}
/* -------------------------------------------------------------------------- */
ModelSolverOptions CouplerSolidContact::getDefaultSolverOptions(
const TimeStepSolverType & type) const {
ModelSolverOptions options;
switch (type) {
case TimeStepSolverType::_dynamic_lumped: {
options.non_linear_solver_type = NonLinearSolverType::_lumped;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_central_difference;
options.solution_type["displacement"] = IntegrationScheme::_acceleration;
break;
}
case TimeStepSolverType::_dynamic: {
options.non_linear_solver_type = NonLinearSolverType::_lumped;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_central_difference;
options.solution_type["displacement"] = IntegrationScheme::_acceleration;
break;
}
case TimeStepSolverType::_static: {
options.non_linear_solver_type =
NonLinearSolverType::_newton_raphson_contact;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_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() {
// computes the internal forces
this->assembleInternalForces();
auto & internal_force = solid->getInternalForce();
auto & external_force = solid->getExternalForce();
auto & contact_force = contact->getInternalForce();
/*auto get_connectivity = [&](auto & slave, auto & master) {
Vector<UInt> master_conn(const_cast<const Mesh &>(mesh).getConnectivity(master));
Vector<UInt> elem_conn(master_conn.size() + 1);
elem_conn[0] = slave;
for (UInt i = 1; i < elem_conn.size(); ++i) {
elem_conn[i] = master_conn[i - 1];
}
return elem_conn;
};
switch(method) {
case _explicit_contact:
case _implicit_contact:
case _explicit_dynamic_contact: {
for (auto & element : contact->getContactElements()) {
for (auto & conn : get_connectivity(element.slave, element.master)) {
for (auto dim : arange(spatial_dimension)) {
external_force(conn, dim) = contact_force(conn, dim);
}
}
}
}
default:
break;
}*/
/* ------------------------------------------------------------------------ */
this->getDOFManager().assembleToResidual("displacement", external_force, 1);
this->getDOFManager().assembleToResidual("displacement", internal_force, 1);
this->getDOFManager().assembleToResidual("displacement", contact_force, 1);
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::assembleResidual(const ID & residual_part) {
AKANTU_DEBUG_IN();
//contact->assembleInternalForces();
auto & internal_force = solid->getInternalForce();
auto & external_force = solid->getExternalForce();
auto & contact_force = contact->getInternalForce();
/*auto get_connectivity = [&](auto & slave, auto & master) {
Vector<UInt> master_conn(const_cast<const Mesh &>(mesh).getConnectivity(master));
Vector<UInt> elem_conn(master_conn.size() + 1);
elem_conn[0] = slave;
for (UInt i = 1; i < elem_conn.size(); ++i) {
elem_conn[i] = master_conn[i - 1];
}
return elem_conn;
};
switch(method) {
case _explicit_dynamic_contact: {
for (auto & element : contact->getContactElements()) {
for (auto & conn : get_connectivity(element.slave, element.master)) {
for (auto dim : arange(spatial_dimension)) {
external_force(conn, dim) = contact_force(conn, dim);
}
}
}
}
default:
break;
}*/
if ("external" == residual_part) {
this->getDOFManager().assembleToResidual("displacement", external_force, 1);
this->getDOFManager().assembleToResidual("displacement", contact_force, 1);
AKANTU_DEBUG_OUT();
return;
}
if ("internal" == residual_part) {
this->getDOFManager().assembleToResidual("displacement", internal_force, 1);
AKANTU_DEBUG_OUT();
return;
}
AKANTU_CUSTOM_EXCEPTION(
debug::SolverCallbackResidualPartUnknown(residual_part));
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::predictor() {
auto & solid_model_solver =
aka::as_type<ModelSolver>(*solid);
solid_model_solver.predictor();
switch (method) {
- case _static:
case _explicit_lumped_mass: {
auto & current_positions = contact->getContactDetector().getPositions();
current_positions.copy(solid->getCurrentPosition());
contact->search();
break;
}
default:
break;
}
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::corrector() {
auto & solid_model_solver =
aka::as_type<ModelSolver>(*solid);
solid_model_solver.corrector();
switch (method) {
case _static:
case _implicit_dynamic: {
auto & current_positions = contact->getContactDetector().getPositions();
current_positions.copy(solid->getCurrentPosition());
contact->search();
break;
}
default:
break;
}
}
/* -------------------------------------------------------------------------- */
MatrixType CouplerSolidContact::getMatrixType(const ID & matrix_id) {
if (matrix_id == "K")
return _symmetric;
if (matrix_id == "M") {
return _symmetric;
}
return _mt_not_defined;
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::assembleMatrix(const ID & matrix_id) {
if (matrix_id == "K") {
this->assembleStiffnessMatrix();
} else if (matrix_id == "M") {
solid->assembleMass();
}
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::assembleLumpedMatrix(const ID & matrix_id) {
if (matrix_id == "M") {
solid->assembleMassLumped();
}
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::beforeSolveStep() {
auto & solid_solver_callback =
aka::as_type<SolverCallback>(*solid);
solid_solver_callback.beforeSolveStep();
auto & contact_solver_callback =
aka::as_type<SolverCallback>(*contact);
contact_solver_callback.beforeSolveStep();
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::afterSolveStep(bool converged) {
auto & solid_solver_callback =
aka::as_type<SolverCallback>(*solid);
solid_solver_callback.afterSolveStep(converged);
auto & contact_solver_callback =
aka::as_type<SolverCallback>(*contact);
contact_solver_callback.afterSolveStep(converged);
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::assembleInternalForces() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the internal forces");
solid->assembleInternalForces();
contact->assembleInternalForces();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::assembleStiffnessMatrix() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the new stiffness matrix");
solid->assembleStiffnessMatrix();
-
+
switch (method) {
case _static:
case _implicit_dynamic: {
- contact->assembleStiffnessMatrix();
+ contact->assembleStiffnessMatrix();
break;
}
default:
break;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::assembleMassLumped() { solid->assembleMassLumped(); }
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::assembleMass() { solid->assembleMass(); }
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::assembleMassLumped(GhostType ghost_type) {
solid->assembleMassLumped(ghost_type);
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::assembleMass(GhostType ghost_type) {
solid->assembleMass(ghost_type);
}
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
/* -------------------------------------------------------------------------- */
std::shared_ptr<dumpers::Field> CouplerSolidContact::createElementalField(
const std::string & field_name, const std::string & group_name,
bool padding_flag, const UInt & spatial_dimension,
const ElementKind & kind) {
std::shared_ptr<dumpers::Field> field;
field = solid->createElementalField(field_name, group_name, padding_flag,
spatial_dimension, kind);
return field;
}
/* -------------------------------------------------------------------------- */
std::shared_ptr<dumpers::Field>
CouplerSolidContact::createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
std::shared_ptr<dumpers::Field> field;
if (field_name == "contact_force" or field_name == "normals" or
field_name == "normal_force" or field_name == "tangential_force" or
field_name == "stick_or_slip" or
field_name == "gaps" or field_name == "previous_gaps" or
field_name == "areas" or field_name == "tangents") {
field = contact->createNodalFieldReal(field_name, group_name, padding_flag);
} else {
field = solid->createNodalFieldReal(field_name, group_name, padding_flag);
}
return field;
}
/* -------------------------------------------------------------------------- */
std::shared_ptr<dumpers::Field>
CouplerSolidContact::createNodalFieldBool(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
std::shared_ptr<dumpers::Field> field;
field = solid->createNodalFieldBool(field_name, group_name, padding_flag);
return field;
}
#else
/* -------------------------------------------------------------------------- */
std::shared_ptr<dumpers::Field>
CouplerSolidContact::createElementalField(const std::string &,
const std::string &, bool,
const UInt &, const ElementKind &) {
return nullptr;
}
/* ----------------------------------------------------------------------- */
std::shared_ptr<dumpers::Field>
CouplerSolidContact::createNodalFieldReal(const std::string &,
const std::string &, bool) {
return nullptr;
}
/*-------------------------------------------------------------------*/
std::shared_ptr<dumpers::Field>
CouplerSolidContact::createNodalFieldBool(const std::string &,
const std::string &, bool) {
return nullptr;
}
#endif
/* --------------------------------------------------------------------------
*/
void CouplerSolidContact::dump(const std::string & dumper_name) {
solid->onDump();
mesh.dump(dumper_name);
}
/* --------------------------------------------------------------------------
*/
void CouplerSolidContact::dump(const std::string & dumper_name, UInt step) {
solid->onDump();
mesh.dump(dumper_name, step);
}
/* -------------------------------------------------------------------------
*/
void CouplerSolidContact::dump(const std::string & dumper_name, Real time,
UInt step) {
solid->onDump();
mesh.dump(dumper_name, time, step);
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::dump() {
solid->onDump();
mesh.dump();
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::dump(UInt step) {
solid->onDump();
mesh.dump(step);
}
/* -------------------------------------------------------------------------- */
void CouplerSolidContact::dump(Real time, UInt step) {
solid->onDump();
mesh.dump(time, step);
}
} // namespace akantu
diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc
index 6bcc5085b..a3040b1f2 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model.cc
@@ -1,1206 +1,1206 @@
/**
* @file solid_mechanics_model.cc
*
* @author Ramin Aghababaei <ramin.aghababaei@epfl.ch>
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Clement Roux <clement.roux@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Tue Jul 27 2010
* @date last modification: Wed Feb 21 2018
*
* @brief Implementation of the SolidMechanicsModel class
*
*
* 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 "solid_mechanics_model.hh"
#include "integrator_gauss.hh"
#include "shape_lagrange.hh"
#include "solid_mechanics_model_tmpl.hh"
#include "communicator.hh"
#include "element_synchronizer.hh"
#include "sparse_matrix.hh"
#include "synchronizer_registry.hh"
#include "dumpable_inline_impl.hh"
#ifdef AKANTU_USE_IOHELPER
#include "dumper_iohelper_paraview.hh"
#endif
#include "material_non_local.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
/**
* A solid mechanics model need a mesh and a dimension to be created. the model
* by it self can not do a lot, the good init functions should be called in
* order to configure the model depending on what we want to do.
*
* @param mesh mesh representing the model we want to simulate
* @param dim spatial dimension of the problem, if dim = 0 (default value) the
* dimension of the problem is assumed to be the on of the mesh
* @param id an id to identify the model
* @param memory_id the id of the memory
* @param model_type this is an internal parameter for inheritance purposes
*/
SolidMechanicsModel::SolidMechanicsModel(
Mesh & mesh, 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),
BoundaryCondition<SolidMechanicsModel>(),
material_index("material index", id, memory_id),
material_local_numbering("material local numbering", id, memory_id) {
AKANTU_DEBUG_IN();
this->registerFEEngineObject<MyFEEngineType>("SolidMechanicsFEEngine", mesh,
Model::spatial_dimension);
#if defined(AKANTU_USE_IOHELPER)
this->mesh.registerDumper<DumperParaview>("solid_mechanics_model", id, true);
this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost,
_ek_regular);
#endif
material_selector = std::make_shared<DefaultMaterialSelector>(material_index);
this->registerDataAccessor(*this);
if (this->mesh.isDistributed()) {
auto & synchronizer = this->mesh.getElementSynchronizer();
this->registerSynchronizer(synchronizer, SynchronizationTag::_material_id);
this->registerSynchronizer(synchronizer, SynchronizationTag::_smm_mass);
this->registerSynchronizer(synchronizer, SynchronizationTag::_smm_stress);
this->registerSynchronizer(synchronizer, SynchronizationTag::_for_dump);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
SolidMechanicsModel::~SolidMechanicsModel() = default;
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::setTimeStep(Real time_step, const ID & solver_id) {
Model::setTimeStep(time_step, solver_id);
#if defined(AKANTU_USE_IOHELPER)
this->mesh.getDumper().setTimeStep(time_step);
#endif
}
/* -------------------------------------------------------------------------- */
/* Initialization */
/* -------------------------------------------------------------------------- */
/**
* This function groups many of the initialization in on function. For most of
* basics case the function should be enough. The functions initialize the
* model, the internal vectors, set them to 0, and depending on the parameters
* it also initialize the explicit or implicit solver.
*
* @param options
* \parblock
* contains the different options to initialize the model
* \li \c analysis_method specify the type of solver to use
* \endparblock
*/
void SolidMechanicsModel::initFullImpl(const ModelOptions & options) {
material_index.initialize(mesh, _element_kind = _ek_not_defined,
_default_value = UInt(-1), _with_nb_element = true);
material_local_numbering.initialize(mesh, _element_kind = _ek_not_defined,
_with_nb_element = true);
Model::initFullImpl(options);
// initialize the materials
if (this->parser.getLastParsedFile() != "") {
this->instantiateMaterials();
this->initMaterials();
}
this->initBC(*this, *displacement, *displacement_increment, *external_force);
}
/* -------------------------------------------------------------------------- */
TimeStepSolverType SolidMechanicsModel::getDefaultSolverType() const {
return TimeStepSolverType::_dynamic_lumped;
}
/* -------------------------------------------------------------------------- */
ModelSolverOptions SolidMechanicsModel::getDefaultSolverOptions(
const TimeStepSolverType & type) const {
ModelSolverOptions options;
switch (type) {
case TimeStepSolverType::_dynamic_lumped: {
options.non_linear_solver_type = NonLinearSolverType::_lumped;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_central_difference;
options.solution_type["displacement"] = IntegrationScheme::_acceleration;
break;
}
case TimeStepSolverType::_static: {
options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_pseudo_time;
options.solution_type["displacement"] = IntegrationScheme::_not_defined;
break;
}
case TimeStepSolverType::_dynamic: {
if (this->method == _explicit_consistent_mass) {
options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_central_difference;
options.solution_type["displacement"] = IntegrationScheme::_acceleration;
} else {
options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_trapezoidal_rule_2;
options.solution_type["displacement"] = IntegrationScheme::_displacement;
}
break;
}
default:
AKANTU_EXCEPTION(type << " is not a valid time step solver type");
}
return options;
}
/* -------------------------------------------------------------------------- */
std::tuple<ID, TimeStepSolverType>
SolidMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) {
switch (method) {
case _explicit_lumped_mass: {
return std::make_tuple("explicit_lumped",
TimeStepSolverType::_dynamic_lumped);
}
case _explicit_consistent_mass: {
return std::make_tuple("explicit", TimeStepSolverType::_dynamic);
}
case _static: {
return std::make_tuple("static", TimeStepSolverType::_static);
}
case _implicit_dynamic: {
return std::make_tuple("implicit", TimeStepSolverType::_dynamic);
}
default:
return std::make_tuple("unknown", TimeStepSolverType::_not_defined);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::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->previous_displacement, spatial_dimension,
"previous_displacement");
this->allocNodalField(this->displacement_increment, spatial_dimension,
"displacement_increment");
this->allocNodalField(this->internal_force, spatial_dimension,
"internal_force");
this->allocNodalField(this->external_force, spatial_dimension,
"external_force");
this->allocNodalField(this->blocked_dofs, spatial_dimension, "blocked_dofs");
this->allocNodalField(this->current_position, spatial_dimension,
"current_position");
// initialize the current positions
this->current_position->copy(this->mesh.getNodes());
/* ------------------------------------------------------------------------ */
if (!dof_manager.hasDOFs("displacement")) {
dof_manager.registerDOFs("displacement", *this->displacement, _dst_nodal);
dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs);
dof_manager.registerDOFsIncrement("displacement",
*this->displacement_increment);
dof_manager.registerDOFsPrevious("displacement",
*this->previous_displacement);
}
/* ------------------------------------------------------------------------ */
// for dynamic
if (time_step_solver_type == TimeStepSolverType::_dynamic ||
time_step_solver_type == TimeStepSolverType::_dynamic_lumped) {
this->allocNodalField(this->velocity, spatial_dimension, "velocity");
this->allocNodalField(this->acceleration, spatial_dimension,
"acceleration");
if (!dof_manager.hasDOFsDerivatives("displacement", 1)) {
dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity);
dof_manager.registerDOFsDerivative("displacement", 2,
*this->acceleration);
}
}
}
/* -------------------------------------------------------------------------- */
/**
* Initialize the model,basically it pre-compute the shapes, shapes derivatives
* and jacobian
*/
void SolidMechanicsModel::initModel() {
/// \todo add the current position as a parameter to initShapeFunctions for
/// large deformation
getFEEngine().initShapeFunctions(_not_ghost);
getFEEngine().initShapeFunctions(_ghost);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleResidual() {
AKANTU_DEBUG_IN();
/* ------------------------------------------------------------------------ */
// computes the internal forces
this->assembleInternalForces();
/* ------------------------------------------------------------------------ */
this->getDOFManager().assembleToResidual("displacement",
*this->external_force, 1);
this->getDOFManager().assembleToResidual("displacement",
*this->internal_force, 1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleResidual(const ID & residual_part) {
AKANTU_DEBUG_IN();
if ("external" == residual_part) {
this->getDOFManager().assembleToResidual("displacement",
*this->external_force, 1);
AKANTU_DEBUG_OUT();
return;
}
if ("internal" == residual_part) {
this->assembleInternalForces();
this->getDOFManager().assembleToResidual("displacement",
*this->internal_force, 1);
AKANTU_DEBUG_OUT();
return;
}
AKANTU_CUSTOM_EXCEPTION(
debug::SolverCallbackResidualPartUnknown(residual_part));
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
MatrixType SolidMechanicsModel::getMatrixType(const ID & matrix_id) {
// \TODO check the materials to know what is the correct answer
if (matrix_id == "C")
return _mt_not_defined;
if (matrix_id == "K") {
auto matrix_type = _unsymmetric;
for (auto & material : materials) {
matrix_type = std::max(matrix_type, material->getMatrixType(matrix_id));
}
}
return _symmetric;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleMatrix(const ID & matrix_id) {
if (matrix_id == "K") {
this->assembleStiffnessMatrix();
} else if (matrix_id == "M") {
this->assembleMass();
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleLumpedMatrix(const ID & matrix_id) {
if (matrix_id == "M") {
this->assembleMassLumped();
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::beforeSolveStep() {
for (auto & material : materials)
material->beforeSolveStep();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::afterSolveStep(bool converged) {
for (auto & material : materials)
material->afterSolveStep(converged);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::predictor() { ++displacement_release; }
/* -------------------------------------------------------------------------- */
-void SolidMechanicsModel::corrector() { ++displacement_release; }
+ void SolidMechanicsModel::corrector() { ++displacement_release; }
/* -------------------------------------------------------------------------- */
/**
* This function computes the internal forces as \f$F_{int} = \int_{\Omega} N
* \sigma d\Omega@\f$
*/
void SolidMechanicsModel::assembleInternalForces() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the internal forces");
this->internal_force->clear();
// compute the stresses of local elements
AKANTU_DEBUG_INFO("Compute local stresses");
for (auto & material : materials) {
material->computeAllStresses(_not_ghost);
}
/* ------------------------------------------------------------------------ */
/* Computation of the non local part */
if (this->non_local_manager)
this->non_local_manager->computeAllNonLocalStresses();
// communicate the stresses
AKANTU_DEBUG_INFO("Send data for residual assembly");
this->asynchronousSynchronize(SynchronizationTag::_smm_stress);
// assemble the forces due to local stresses
AKANTU_DEBUG_INFO("Assemble residual for local elements");
for (auto & material : materials) {
material->assembleInternalForces(_not_ghost);
}
// finalize communications
AKANTU_DEBUG_INFO("Wait distant stresses");
this->waitEndSynchronize(SynchronizationTag::_smm_stress);
// assemble the stresses due to ghost elements
AKANTU_DEBUG_INFO("Assemble residual for ghost elements");
for (auto & material : materials) {
material->assembleInternalForces(_ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleStiffnessMatrix() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the new stiffness matrix.");
// Check if materials need to recompute the matrix
bool need_to_reassemble = true;
/*bool need_to_reassemble = false;
for (auto & material : materials) {
need_to_reassemble |= material->hasMatrixChanged("K");
}*/
if (need_to_reassemble) {
this->getDOFManager().getMatrix("K").clear();
// call compute stiffness matrix on each local elements
for (auto & material : materials) {
material->assembleStiffnessMatrix(_not_ghost);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateCurrentPosition() {
if (this->current_position_release == this->displacement_release)
return;
this->current_position->copy(this->mesh.getNodes());
auto cpos_it = this->current_position->begin(Model::spatial_dimension);
auto cpos_end = this->current_position->end(Model::spatial_dimension);
auto disp_it = this->displacement->begin(Model::spatial_dimension);
for (; cpos_it != cpos_end; ++cpos_it, ++disp_it) {
*cpos_it += *disp_it;
}
this->current_position_release = this->displacement_release;
}
/* -------------------------------------------------------------------------- */
const Array<Real> & SolidMechanicsModel::getCurrentPosition() {
this->updateCurrentPosition();
return *this->current_position;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateDataForNonLocalCriterion(
ElementTypeMapReal & criterion) {
const ID field_name = criterion.getName();
for (auto & material : materials) {
if (!material->isInternal<Real>(field_name, _ek_regular))
continue;
for (auto ghost_type : ghost_types) {
material->flattenInternal(field_name, criterion, ghost_type, _ek_regular);
}
}
}
/* -------------------------------------------------------------------------- */
/* Information */
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getStableTimeStep() {
AKANTU_DEBUG_IN();
Real min_dt = getStableTimeStep(_not_ghost);
/// reduction min over all processors
mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min);
AKANTU_DEBUG_OUT();
return min_dt;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) {
AKANTU_DEBUG_IN();
Real min_dt = std::numeric_limits<Real>::max();
this->updateCurrentPosition();
Element elem;
elem.ghost_type = ghost_type;
for (auto type :
mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_regular)) {
elem.type = type;
UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type);
UInt nb_element = mesh.getNbElement(type);
auto mat_indexes = material_index(type, ghost_type).begin();
auto mat_loc_num = material_local_numbering(type, ghost_type).begin();
Array<Real> X(0, nb_nodes_per_element * Model::spatial_dimension);
FEEngine::extractNodalToElementField(mesh, *current_position, X, type,
_not_ghost);
auto X_el = X.begin(Model::spatial_dimension, nb_nodes_per_element);
for (UInt el = 0; el < nb_element;
++el, ++X_el, ++mat_indexes, ++mat_loc_num) {
elem.element = *mat_loc_num;
Real el_h = getFEEngine().getElementInradius(*X_el, type);
Real el_c = this->materials[*mat_indexes]->getCelerity(elem);
Real el_dt = el_h / el_c;
min_dt = std::min(min_dt, el_dt);
}
}
AKANTU_DEBUG_OUT();
return min_dt;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getKineticEnergy() {
AKANTU_DEBUG_IN();
Real ekin = 0.;
UInt nb_nodes = mesh.getNbNodes();
if (this->getDOFManager().hasLumpedMatrix("M")) {
auto m_it = this->mass->begin(Model::spatial_dimension);
auto m_end = this->mass->end(Model::spatial_dimension);
auto v_it = this->velocity->begin(Model::spatial_dimension);
for (UInt n = 0; m_it != m_end; ++n, ++m_it, ++v_it) {
const auto & v = *v_it;
const auto & m = *m_it;
Real mv2 = 0.;
auto is_local_node = mesh.isLocalOrMasterNode(n);
// bool is_not_pbc_slave_node = !isPBCSlaveNode(n);
auto count_node = is_local_node; // && is_not_pbc_slave_node;
if (count_node) {
for (UInt i = 0; i < Model::spatial_dimension; ++i) {
if (m(i) > std::numeric_limits<Real>::epsilon())
mv2 += v(i) * v(i) * m(i);
}
}
ekin += mv2;
}
} else if (this->getDOFManager().hasMatrix("M")) {
Array<Real> Mv(nb_nodes, Model::spatial_dimension);
this->getDOFManager().assembleMatMulVectToArray("displacement", "M",
*this->velocity, Mv);
for (auto && data : zip(arange(nb_nodes), make_view(Mv, spatial_dimension),
make_view(*this->velocity, spatial_dimension))) {
ekin += std::get<2>(data).dot(std::get<1>(data)) *
mesh.isLocalOrMasterNode(std::get<0>(data));
}
} else {
AKANTU_ERROR("No function called to assemble the mass matrix.");
}
mesh.getCommunicator().allReduce(ekin, SynchronizerOperation::_sum);
AKANTU_DEBUG_OUT();
return ekin * .5;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getKineticEnergy(const ElementType & type,
UInt index) {
AKANTU_DEBUG_IN();
UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type);
Array<Real> vel_on_quad(nb_quadrature_points, Model::spatial_dimension);
Array<UInt> filter_element(1, 1, index);
getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad,
Model::spatial_dimension, type,
_not_ghost, filter_element);
auto vit = vel_on_quad.begin(Model::spatial_dimension);
auto vend = vel_on_quad.end(Model::spatial_dimension);
Vector<Real> rho_v2(nb_quadrature_points);
Real rho = materials[material_index(type)(index)]->getRho();
for (UInt q = 0; vit != vend; ++vit, ++q) {
rho_v2(q) = rho * vit->dot(*vit);
}
AKANTU_DEBUG_OUT();
return .5 * getFEEngine().integrate(rho_v2, type, index);
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getExternalWork() {
AKANTU_DEBUG_IN();
auto ext_force_it = external_force->begin(Model::spatial_dimension);
auto int_force_it = internal_force->begin(Model::spatial_dimension);
auto boun_it = blocked_dofs->begin(Model::spatial_dimension);
decltype(ext_force_it) incr_or_velo_it;
if (this->method == _static) {
incr_or_velo_it =
this->displacement_increment->begin(Model::spatial_dimension);
} else {
incr_or_velo_it = this->velocity->begin(Model::spatial_dimension);
}
Real work = 0.;
UInt nb_nodes = this->mesh.getNbNodes();
for (UInt n = 0; n < nb_nodes;
++n, ++ext_force_it, ++int_force_it, ++boun_it, ++incr_or_velo_it) {
const auto & int_force = *int_force_it;
const auto & ext_force = *ext_force_it;
const auto & boun = *boun_it;
const auto & incr_or_velo = *incr_or_velo_it;
bool is_local_node = this->mesh.isLocalOrMasterNode(n);
// bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n);
bool count_node = is_local_node; // && is_not_pbc_slave_node;
if (count_node) {
for (UInt i = 0; i < Model::spatial_dimension; ++i) {
if (boun(i))
work -= int_force(i) * incr_or_velo(i);
else
work += ext_force(i) * incr_or_velo(i);
}
}
}
mesh.getCommunicator().allReduce(work, SynchronizerOperation::_sum);
if (this->method != _static)
work *= this->getTimeStep();
AKANTU_DEBUG_OUT();
return work;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getEnergy(const std::string & energy_id) {
AKANTU_DEBUG_IN();
if (energy_id == "kinetic") {
return getKineticEnergy();
} else if (energy_id == "external work") {
return getExternalWork();
}
Real energy = 0.;
for (auto & material : materials)
energy += material->getEnergy(energy_id);
/// reduction sum over all processors
mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum);
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getEnergy(const std::string & energy_id,
const ElementType & type, UInt index) {
AKANTU_DEBUG_IN();
if (energy_id == "kinetic") {
return getKineticEnergy(type, index);
}
UInt mat_index = this->material_index(type, _not_ghost)(index);
UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index);
Real energy =
this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num);
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onElementsAdded(const Array<Element> & element_list,
const NewElementsEvent & event) {
AKANTU_DEBUG_IN();
this->material_index.initialize(mesh, _element_kind = _ek_not_defined,
_with_nb_element = true,
_default_value = UInt(-1));
this->material_local_numbering.initialize(
mesh, _element_kind = _ek_not_defined, _with_nb_element = true,
_default_value = UInt(-1));
ElementTypeMapArray<UInt> filter("new_element_filter", this->getID(),
this->getMemoryID());
for (auto & elem : element_list) {
if (mesh.getSpatialDimension(elem.type) != spatial_dimension)
continue;
if (!filter.exists(elem.type, elem.ghost_type))
filter.alloc(0, 1, elem.type, elem.ghost_type);
filter(elem.type, elem.ghost_type).push_back(elem.element);
}
// this fails in parallel if the event is sent on facet between constructor
// and initFull \todo: to debug...
this->assignMaterialToElements(&filter);
for (auto & material : materials)
material->onElementsAdded(element_list, event);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onElementsRemoved(
const Array<Element> & element_list,
const ElementTypeMapArray<UInt> & new_numbering,
const RemovedElementsEvent & event) {
for (auto & material : materials) {
material->onElementsRemoved(element_list, new_numbering, event);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onNodesAdded(const Array<UInt> & nodes_list,
const NewNodesEvent & event) {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
if (displacement) {
displacement->resize(nb_nodes, 0.);
++displacement_release;
}
if (mass)
mass->resize(nb_nodes, 0.);
if (velocity)
velocity->resize(nb_nodes, 0.);
if (acceleration)
acceleration->resize(nb_nodes, 0.);
if (external_force)
external_force->resize(nb_nodes, 0.);
if (internal_force)
internal_force->resize(nb_nodes, 0.);
if (blocked_dofs)
blocked_dofs->resize(nb_nodes, 0.);
if (current_position)
current_position->resize(nb_nodes, 0.);
if (previous_displacement)
previous_displacement->resize(nb_nodes, 0.);
if (displacement_increment)
displacement_increment->resize(nb_nodes, 0.);
for (auto & material : materials) {
material->onNodesAdded(nodes_list, event);
}
need_to_reassemble_lumped_mass = true;
need_to_reassemble_mass = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onNodesRemoved(const Array<UInt> & /*element_list*/,
const Array<UInt> & new_numbering,
const RemovedNodesEvent & /*event*/) {
if (displacement) {
mesh.removeNodesFromArray(*displacement, new_numbering);
++displacement_release;
}
if (mass)
mesh.removeNodesFromArray(*mass, new_numbering);
if (velocity)
mesh.removeNodesFromArray(*velocity, new_numbering);
if (acceleration)
mesh.removeNodesFromArray(*acceleration, new_numbering);
if (internal_force)
mesh.removeNodesFromArray(*internal_force, new_numbering);
if (external_force)
mesh.removeNodesFromArray(*external_force, new_numbering);
if (blocked_dofs)
mesh.removeNodesFromArray(*blocked_dofs, new_numbering);
// if (increment_acceleration)
// mesh.removeNodesFromArray(*increment_acceleration, new_numbering);
if (displacement_increment)
mesh.removeNodesFromArray(*displacement_increment, new_numbering);
if (previous_displacement)
mesh.removeNodesFromArray(*previous_displacement, new_numbering);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::printself(std::ostream & stream, int indent) const {
std::string space(indent, AKANTU_INDENT);
stream << space << "Solid 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 << " ]" << std::endl;
stream << space << " + nodals information [" << std::endl;
displacement->printself(stream, indent + 2);
if (velocity)
velocity->printself(stream, indent + 2);
if (acceleration)
acceleration->printself(stream, indent + 2);
if (mass)
mass->printself(stream, indent + 2);
external_force->printself(stream, indent + 2);
internal_force->printself(stream, indent + 2);
blocked_dofs->printself(stream, indent + 2);
stream << space << " ]" << std::endl;
stream << space << " + material information [" << std::endl;
material_index.printself(stream, indent + 2);
stream << space << " ]" << std::endl;
stream << space << " + materials [" << std::endl;
for (auto & material : materials)
material->printself(stream, indent + 2);
stream << space << " ]" << std::endl;
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initializeNonLocal() {
this->non_local_manager->synchronize(*this, SynchronizationTag::_material_id);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::insertIntegrationPointsInNeighborhoods(
const GhostType & ghost_type) {
for (auto & mat : materials) {
MaterialNonLocalInterface * mat_non_local;
if ((mat_non_local =
dynamic_cast<MaterialNonLocalInterface *>(mat.get())) == nullptr)
continue;
ElementTypeMapArray<Real> quadrature_points_coordinates(
"quadrature_points_coordinates_tmp_nl", this->id, this->memory_id);
quadrature_points_coordinates.initialize(this->getFEEngine(),
_nb_component = spatial_dimension,
_ghost_type = ghost_type);
for (auto & type : quadrature_points_coordinates.elementTypes(
Model::spatial_dimension, ghost_type)) {
this->getFEEngine().computeIntegrationPointsCoordinates(
quadrature_points_coordinates(type, ghost_type), type, ghost_type);
}
mat_non_local->initMaterialNonLocal();
mat_non_local->insertIntegrationPointsInNeighborhoods(
ghost_type, quadrature_points_coordinates);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::computeNonLocalStresses(
const GhostType & ghost_type) {
for (auto & mat : materials) {
if (not aka::is_of_type<MaterialNonLocalInterface>(*mat))
continue;
auto & mat_non_local = dynamic_cast<MaterialNonLocalInterface &>(*mat);
mat_non_local.computeNonLocalStresses(ghost_type);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateLocalInternal(
ElementTypeMapReal & internal_flat, const GhostType & ghost_type,
const ElementKind & kind) {
const ID field_name = internal_flat.getName();
for (auto & material : materials) {
if (material->isInternal<Real>(field_name, kind))
material->flattenInternal(field_name, internal_flat, ghost_type, kind);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateNonLocalInternal(
ElementTypeMapReal & internal_flat, const GhostType & ghost_type,
const ElementKind & kind) {
const ID field_name = internal_flat.getName();
for (auto & mat : materials) {
if (not aka::is_of_type<MaterialNonLocalInterface>(*mat))
continue;
auto & mat_non_local = dynamic_cast<MaterialNonLocalInterface &>(*mat);
mat_non_local.updateNonLocalInternals(internal_flat, field_name, ghost_type,
kind);
}
}
/* -------------------------------------------------------------------------- */
FEEngine & SolidMechanicsModel::getFEEngineBoundary(const ID & name) {
return getFEEngineClassBoundary<MyFEEngineType>(name);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::splitElementByMaterial(
const Array<Element> & elements,
std::vector<Array<Element>> & elements_per_mat) const {
for (const auto & el : elements) {
Element mat_el = el;
mat_el.element = this->material_local_numbering(el);
elements_per_mat[this->material_index(el)].push_back(mat_el);
}
}
/* -------------------------------------------------------------------------- */
UInt SolidMechanicsModel::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);
}
switch (tag) {
case SynchronizationTag::_material_id: {
size += elements.size() * sizeof(UInt);
break;
}
case SynchronizationTag::_smm_mass: {
size += nb_nodes_per_element * sizeof(Real) *
Model::spatial_dimension; // mass vector
break;
}
case SynchronizationTag::_smm_for_gradu: {
size += nb_nodes_per_element * Model::spatial_dimension *
sizeof(Real); // displacement
break;
}
case SynchronizationTag::_smm_boundary: {
// force, displacement, boundary
size += nb_nodes_per_element * Model::spatial_dimension *
(2 * sizeof(Real) + sizeof(bool));
break;
}
case SynchronizationTag::_for_dump: {
// displacement, velocity, acceleration, residual, force
size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real) * 5;
break;
}
default: {
}
}
if (tag != SynchronizationTag::_material_id) {
splitByMaterial(elements, [&](auto && mat, auto && elements) {
size += mat.getNbData(elements, tag);
});
}
AKANTU_DEBUG_OUT();
return size;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::packData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) const {
AKANTU_DEBUG_IN();
switch (tag) {
case SynchronizationTag::_material_id: {
this->packElementalDataHelper(material_index, buffer, elements, false,
getFEEngine());
break;
}
case SynchronizationTag::_smm_mass: {
packNodalDataHelper(*mass, buffer, elements, mesh);
break;
}
case SynchronizationTag::_smm_for_gradu: {
packNodalDataHelper(*displacement, buffer, elements, mesh);
break;
}
case SynchronizationTag::_for_dump: {
packNodalDataHelper(*displacement, buffer, elements, mesh);
packNodalDataHelper(*velocity, buffer, elements, mesh);
packNodalDataHelper(*acceleration, buffer, elements, mesh);
packNodalDataHelper(*internal_force, buffer, elements, mesh);
packNodalDataHelper(*external_force, buffer, elements, mesh);
break;
}
case SynchronizationTag::_smm_boundary: {
packNodalDataHelper(*external_force, buffer, elements, mesh);
packNodalDataHelper(*velocity, buffer, elements, mesh);
packNodalDataHelper(*blocked_dofs, buffer, elements, mesh);
break;
}
default: {
}
}
if (tag != SynchronizationTag::_material_id) {
splitByMaterial(elements, [&](auto && mat, auto && elements) {
mat.packData(buffer, elements, tag);
});
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) {
AKANTU_DEBUG_IN();
switch (tag) {
case SynchronizationTag::_material_id: {
for (auto && element : elements) {
UInt recv_mat_index;
buffer >> recv_mat_index;
UInt & mat_index = material_index(element);
if (mat_index != UInt(-1))
continue;
// add ghosts element to the correct material
mat_index = recv_mat_index;
UInt index = materials[mat_index]->addElement(element);
material_local_numbering(element) = index;
}
break;
}
case SynchronizationTag::_smm_mass: {
unpackNodalDataHelper(*mass, buffer, elements, mesh);
break;
}
case SynchronizationTag::_smm_for_gradu: {
unpackNodalDataHelper(*displacement, buffer, elements, mesh);
break;
}
case SynchronizationTag::_for_dump: {
unpackNodalDataHelper(*displacement, buffer, elements, mesh);
unpackNodalDataHelper(*velocity, buffer, elements, mesh);
unpackNodalDataHelper(*acceleration, buffer, elements, mesh);
unpackNodalDataHelper(*internal_force, buffer, elements, mesh);
unpackNodalDataHelper(*external_force, buffer, elements, mesh);
break;
}
case SynchronizationTag::_smm_boundary: {
unpackNodalDataHelper(*external_force, buffer, elements, mesh);
unpackNodalDataHelper(*velocity, buffer, elements, mesh);
unpackNodalDataHelper(*blocked_dofs, buffer, elements, mesh);
break;
}
default: {
}
}
if (tag != SynchronizationTag::_material_id) {
splitByMaterial(elements, [&](auto && mat, auto && elements) {
mat.unpackData(buffer, elements, tag);
});
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
UInt SolidMechanicsModel::getNbData(const Array<UInt> & dofs,
const SynchronizationTag & tag) const {
AKANTU_DEBUG_IN();
UInt size = 0;
// UInt nb_nodes = mesh.getNbNodes();
switch (tag) {
case SynchronizationTag::_smm_uv: {
size += sizeof(Real) * Model::spatial_dimension * 2;
break;
}
case SynchronizationTag::_smm_res: {
size += sizeof(Real) * Model::spatial_dimension;
break;
}
case SynchronizationTag::_smm_mass: {
size += sizeof(Real) * Model::spatial_dimension;
break;
}
case SynchronizationTag::_for_dump: {
size += sizeof(Real) * Model::spatial_dimension * 5;
break;
}
default: {
AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
AKANTU_DEBUG_OUT();
return size * dofs.size();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::packData(CommunicationBuffer & buffer,
const Array<UInt> & dofs,
const SynchronizationTag & tag) const {
AKANTU_DEBUG_IN();
switch (tag) {
case SynchronizationTag::_smm_uv: {
packDOFDataHelper(*displacement, buffer, dofs);
packDOFDataHelper(*velocity, buffer, dofs);
break;
}
case SynchronizationTag::_smm_res: {
packDOFDataHelper(*internal_force, buffer, dofs);
break;
}
case SynchronizationTag::_smm_mass: {
packDOFDataHelper(*mass, buffer, dofs);
break;
}
case SynchronizationTag::_for_dump: {
packDOFDataHelper(*displacement, buffer, dofs);
packDOFDataHelper(*velocity, buffer, dofs);
packDOFDataHelper(*acceleration, buffer, dofs);
packDOFDataHelper(*internal_force, buffer, dofs);
packDOFDataHelper(*external_force, buffer, dofs);
break;
}
default: {
AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer,
const Array<UInt> & dofs,
const SynchronizationTag & tag) {
AKANTU_DEBUG_IN();
switch (tag) {
case SynchronizationTag::_smm_uv: {
unpackDOFDataHelper(*displacement, buffer, dofs);
unpackDOFDataHelper(*velocity, buffer, dofs);
break;
}
case SynchronizationTag::_smm_res: {
unpackDOFDataHelper(*internal_force, buffer, dofs);
break;
}
case SynchronizationTag::_smm_mass: {
unpackDOFDataHelper(*mass, buffer, dofs);
break;
}
case SynchronizationTag::_for_dump: {
unpackDOFDataHelper(*displacement, buffer, dofs);
unpackDOFDataHelper(*velocity, buffer, dofs);
unpackDOFDataHelper(*acceleration, buffer, dofs);
unpackDOFDataHelper(*internal_force, buffer, dofs);
unpackDOFDataHelper(*external_force, buffer, dofs);
break;
}
default: {
AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
AKANTU_DEBUG_OUT();
}
} // namespace akantu