diff --git a/extra_packages/extra-materials b/extra_packages/extra-materials
index 27835edcf..b6e70a3df 160000
--- a/extra_packages/extra-materials
+++ b/extra_packages/extra-materials
@@ -1 +1 @@
-Subproject commit 27835edcf205dedb552ddb9d82d3bc9e39b0156a
+Subproject commit b6e70a3dfb54b4021a10b37a9e4e6c26a64f9854
diff --git a/extra_packages/parallel-cohesive-element b/extra_packages/parallel-cohesive-element
index 93c04c339..c46ced178 160000
--- a/extra_packages/parallel-cohesive-element
+++ b/extra_packages/parallel-cohesive-element
@@ -1 +1 @@
-Subproject commit 93c04c33905f63fc4addef7160fa400409b154fb
+Subproject commit c46ced1786aec62d4e082dd3624eb539a96413f4
diff --git a/python/swig/aka_common.i b/python/swig/aka_common.i
index 6079d8ef9..fd9c3b4ec 100644
--- a/python/swig/aka_common.i
+++ b/python/swig/aka_common.i
@@ -1,25 +1,30 @@
%{
#include "aka_common.hh"
%}
namespace akantu {
%ignore getStaticParser;
%ignore getUserParser;
%ignore initialize(int & argc, char ** & argv);
%ignore initialize(const std::string & input_file, int & argc, char ** & argv);
}
%inline %{
namespace akantu {
void initialize(const std::string & input_file) {
int argc = 0;
char ** argv = NULL;
initialize(input_file, argc, argv);
}
+ void initialize() {
+ int argc = 0;
+ char ** argv = NULL;
+ initialize(argc, argv);
+ }
}
%}
%include "aka_common.hh"
diff --git a/python/swig/akantu.i b/python/swig/akantu.i
index 31ab927fe..d9232e07c 100644
--- a/python/swig/akantu.i
+++ b/python/swig/akantu.i
@@ -1,26 +1,29 @@
%module akantu
%include "stl.i"
#define __attribute__(x)
%ignore akantu::operator <<;
%include "aka_common.i"
%include "aka_array.i"
%define print_self(MY_CLASS)
%extend akantu::MY_CLASS {
std::string __str__() {
std::stringstream sstr;
sstr << *($self);
return sstr.str();
}
}
%enddef
%include "mesh.i"
%include "mesh_utils.i"
%include "model.i"
%include "solid_mechanics_model.i"
-
+#if defined(AKANTU_STRUCTURAL_MECHANICS)
+%include "load_functions.i"
+%include "structural_mechanics_model.i"
+#endif
diff --git a/python/swig/load_functions.i b/python/swig/load_functions.i
new file mode 100644
index 000000000..046c3bfe0
--- /dev/null
+++ b/python/swig/load_functions.i
@@ -0,0 +1,14 @@
+%inline %{
+ namespace akantu {
+
+ static void lin_load(double * position, double * load,
+ __attribute__ ((unused)) Real * normal,
+ __attribute__ ((unused)) UInt surface_id) {
+
+ memset(load,0,sizeof(Real)*3);
+ if (position[0]<=10){
+ load[1]= -6000;
+ }
+ }
+ }
+ %}
diff --git a/python/swig/mesh.i b/python/swig/mesh.i
index 283f68965..d5740a174 100644
--- a/python/swig/mesh.i
+++ b/python/swig/mesh.i
@@ -1,95 +1,113 @@
%{
#include "mesh.hh"
#include "node_group.hh"
#include "solid_mechanics_model.hh"
using akantu::Vector;
using akantu::ElementTypeMapArray;
using akantu::MatrixProxy;
using akantu::Matrix;
using akantu::UInt;
using akantu::Real;
using akantu::Array;
using akantu::SolidMechanicsModel;
%}
namespace akantu {
%ignore NewNodesEvent;
%ignore RemovedNodesEvent;
%ignore NewElementsEvent;
%ignore RemovedElementsEvent;
%ignore MeshEventHandler;
%ignore MeshEvent< UInt >;
%ignore MeshEvent< Element >;
%ignore Mesh::extractNodalCoordinatesFromPBCElement;
%ignore Mesh::getGroupDumer;
}
print_self(Mesh)
+
+%extend akantu::Mesh {
+
+ void resizeMesh(UInt nb_nodes, UInt nb_element, const ElementType & type) {
+
+ Array<Real> & nodes = const_cast<Array<Real> &>($self->getNodes());
+ nodes.resize(nb_nodes);
+
+ $self->addConnectivityType(type);
+ Array<UInt> & connectivity = const_cast<Array<UInt> &>($self->getConnectivity(type));
+ connectivity.resize(nb_element);
+
+
+
+ }
+
+}
+
%extend akantu::GroupManager {
void createGroupsFromStringMeshData(const std::string & dataset_name) {
$self->createGroupsFromMeshData<std::string>(dataset_name);
}
void createGroupsFromUIntMeshData(const std::string & dataset_name) {
$self->createGroupsFromMeshData<akantu::UInt>(dataset_name);
}
}
%extend akantu::NodeGroup {
akantu::Array<akantu::Real> & getGroupedNodes(akantu::Array<akantu::Real, true> & surface_array, Mesh & mesh) {
akantu::Array<akantu::UInt> group_node = $self->getNodes();
akantu::Array<akantu::Real> & full_array = mesh.getNodes();
surface_array.resize(group_node.getSize());
for (UInt i = 0; i < group_node.getSize(); ++i) {
for (UInt cmp = 0; cmp < full_array.getNbComponent(); ++cmp) {
surface_array(i,cmp) = full_array(group_node(i),cmp);
}
}
akantu::Array<akantu::Real> & res(surface_array);
return res;
}
akantu::Array<akantu::Real> & getGroupedArray(akantu::Array<akantu::Real, true> & surface_array, akantu::SolidMechanicsModel & model, int type) {
akantu::Array<akantu::Real> * full_array;
switch (type) {
case 0 : full_array = new akantu::Array<akantu::Real>(model.getDisplacement());
break;
case 1 : full_array = new akantu::Array<akantu::Real>(model.getVelocity());
break;
case 2 : full_array = new akantu::Array<akantu::Real>(model.getForce());
break;
}
akantu::Array<akantu::UInt> group_node = $self->getNodes();
surface_array.resize(group_node.getSize());
for (UInt i = 0; i < group_node.getSize(); ++i) {
for (UInt cmp = 0; cmp < full_array->getNbComponent(); ++cmp) {
surface_array(i,cmp) = (*full_array)(group_node(i),cmp);
}
}
akantu::Array<akantu::Real> & res(surface_array);
return res;
}
}
%include "group_manager.hh"
%include "element_group.hh"
%include "node_group.hh"
%include "mesh.hh"
diff --git a/python/swig/structural_mechanics_model.i b/python/swig/structural_mechanics_model.i
new file mode 100644
index 000000000..c6f45e515
--- /dev/null
+++ b/python/swig/structural_mechanics_model.i
@@ -0,0 +1,40 @@
+%{
+ #include "structural_mechanics_model.hh"
+ #include "sparse_matrix.hh"
+%}
+
+namespace akantu {
+ %ignore StructuralMechanicsModel::initArrays;
+ %ignore StructuralMechanicsModel::initModel;
+
+ %ignore StructuralMechanicsModel::initSolver;
+ %ignore StructuralMechanicsModel::initImplicit;
+
+ static void lin_load(double * position, double * load,
+ __attribute__ ((unused)) Real * normal,
+ __attribute__ ((unused)) UInt surface_id) {
+
+ memset(load,0,sizeof(Real)*3);
+ if (position[0]<=10){
+ load[1]= -6000;
+ }
+ }
+}
+
+%include "dumpable.hh"
+%include "structural_mechanics_model.hh"
+
+%extend akantu::StructuralMechanicsModel {
+
+ bool solveStep(Real tolerance, UInt max_iteration) {
+
+ return $self->solveStep<akantu::SolveConvergenceMethod::_scm_newton_raphson_tangent, akantu::SolveConvergenceCriteria::_scc_residual>(tolerance, max_iteration);
+
+ }
+
+ void computeForcesFromFunctionBeam2d(BoundaryFunctionType function_type) {
+
+ $self->computeForcesFromFunction<akantu::ElementType::_bernoulli_beam_2>(akantu::lin_load, function_type);
+
+ }
+ }
diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc
index 1cbf916c7..352ffd4b6 100644
--- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc
+++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc
@@ -1,655 +1,656 @@
/**
* @file material_cohesive_linear.cc
*
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Tue May 08 2012
* @date last modification: Thu Aug 07 2014
*
* @brief Linear irreversible cohesive law of mixed mode loading with
* random stress definition for extrinsic type
*
* @section LICENSE
*
* Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include <numeric>
/* -------------------------------------------------------------------------- */
#include "material_cohesive_linear.hh"
#include "solid_mechanics_model_cohesive.hh"
#include "sparse_matrix.hh"
#include "dof_synchronizer.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
MaterialCohesiveLinear<spatial_dimension>::MaterialCohesiveLinear(SolidMechanicsModel & model,
const ID & id) :
MaterialCohesive(model,id),
sigma_c_eff("sigma_c_eff", *this),
delta_c_eff("delta_c_eff", *this),
insertion_stress("insertion_stress", *this) {
AKANTU_DEBUG_IN();
this->registerParam("beta" , beta , 0. ,
_pat_parsable | _pat_readable,
"Beta parameter" );
this->registerParam("G_c" , G_c , 0. ,
_pat_parsable | _pat_readable,
"Mode I fracture energy" );
this->registerParam("penalty", penalty, 0. ,
_pat_parsable | _pat_readable,
"Penalty coefficient" );
this->registerParam("volume_s", volume_s, 0. ,
_pat_parsable | _pat_readable,
"Reference volume for sigma_c scaling");
this->registerParam("m_s", m_s, 1. ,
_pat_parsable | _pat_readable,
"Weibull exponent for sigma_c scaling");
this->registerParam("kappa" , kappa , 1. ,
_pat_parsable | _pat_readable,
"Kappa parameter");
// if (!model->isExplicit())
use_previous_delta_max = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialCohesiveLinear<spatial_dimension>::initMaterial() {
AKANTU_DEBUG_IN();
MaterialCohesive::initMaterial();
/// compute scalars
beta2_kappa2 = beta * beta/kappa/kappa;
beta2_kappa = beta * beta/kappa;
if (Math::are_float_equal(beta, 0))
beta2_inv = 0;
else
beta2_inv = 1./beta/beta;
sigma_c_eff.initialize(1);
delta_c_eff.initialize(1);
insertion_stress.initialize(spatial_dimension);
if (!Math::are_float_equal(delta_c, 0.))
delta_c_eff.setDefaultValue(delta_c);
else
delta_c_eff.setDefaultValue(2 * G_c / sigma_c);
if (model->getIsExtrinsic()) scaleInsertionTraction();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialCohesiveLinear<spatial_dimension>::scaleInsertionTraction() {
AKANTU_DEBUG_IN();
// do nothing if volume_s hasn't been specified by the user
if (Math::are_float_equal(volume_s, 0.)) return;
const Mesh & mesh_facets = model->getMeshFacets();
const FEEngine & fe_engine = model->getFEEngine();
const FEEngine & fe_engine_facet = model->getFEEngine("FacetsFEEngine");
// loop over facet type
Mesh::type_iterator first = mesh_facets.firstType(spatial_dimension - 1);
Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1);
Real base_sigma_c = sigma_c;
for(;first != last; ++first) {
ElementType type_facet = *first;
const Array< std::vector<Element> > & facet_to_element
= mesh_facets.getElementToSubelement(type_facet);
UInt nb_facet = facet_to_element.getSize();
UInt nb_quad_per_facet = fe_engine_facet.getNbQuadraturePoints(type_facet);
// iterator to modify sigma_c for all the quadrature points of a facet
Array<Real>::vector_iterator sigma_c_iterator
= sigma_c(type_facet).begin_reinterpret(nb_quad_per_facet, nb_facet);
for (UInt f = 0; f < nb_facet; ++f, ++sigma_c_iterator) {
const std::vector<Element> & element_list = facet_to_element(f);
// compute bounding volume
Real volume = 0;
std::vector<Element>::const_iterator elem = element_list.begin();
std::vector<Element>::const_iterator elem_end = element_list.end();
for (; elem != elem_end; ++elem) {
if (*elem == ElementNull) continue;
// unit vector for integration in order to obtain the volume
UInt nb_quadrature_points = fe_engine.getNbQuadraturePoints(elem->type);
Vector<Real> unit_vector(nb_quadrature_points, 1);
volume += fe_engine.integrate(unit_vector, elem->type,
elem->element, elem->ghost_type);
}
// scale sigma_c
*sigma_c_iterator -= base_sigma_c;
*sigma_c_iterator *= std::pow(volume_s / volume, 1. / m_s);
*sigma_c_iterator += base_sigma_c;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialCohesiveLinear<spatial_dimension>::checkInsertion() {
AKANTU_DEBUG_IN();
const Mesh & mesh_facets = model->getMeshFacets();
CohesiveElementInserter & inserter = model->getElementInserter();
Real tolerance = Math::getTolerance();
Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1);
Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1);
for (; it != last; ++it) {
ElementType type_facet = *it;
ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet);
const Array<bool> & facets_check = inserter.getCheckFacets(type_facet);
Array<bool> & f_insertion = inserter.getInsertionFacets(type_facet);
Array<UInt> & f_filter = facet_filter(type_facet);
Array<Real> & sig_c_eff = sigma_c_eff(type_cohesive);
Array<Real> & del_c = delta_c_eff(type_cohesive);
Array<Real> & ins_stress = insertion_stress(type_cohesive);
Array<Real> & trac_old = tractions_old(type_cohesive);
const Array<Real> & f_stress = model->getStressOnFacets(type_facet);
const Array<Real> & sigma_lim = sigma_c(type_facet);
Real max_ratio = 0.;
UInt index_f = 0;
UInt index_filter = 0;
UInt nn = 0;
//Real sigma_c_tmp = 0.0;
UInt nb_quad_facet = model->getFEEngine("FacetsFEEngine").getNbQuadraturePoints(type_facet);
UInt nb_facet = f_filter.getSize();
if (nb_facet == 0) continue;
Array<Real>::const_iterator<Real> sigma_lim_it = sigma_lim.begin();
Matrix<Real> stress_tmp(spatial_dimension, spatial_dimension);
Matrix<Real> normal_traction(spatial_dimension, nb_quad_facet);
Vector<Real> stress_check(nb_quad_facet);
UInt sp2 = spatial_dimension * spatial_dimension;
const Array<Real> & tangents = model->getTangents(type_facet);
const Array<Real> & normals
= model->getFEEngine("FacetsFEEngine").getNormalsOnQuadPoints(type_facet);
Array<Real>::const_vector_iterator normal_begin = normals.begin(spatial_dimension);
Array<Real>::const_vector_iterator tangent_begin = tangents.begin(tangents.getNbComponent());
Array<Real>::const_matrix_iterator facet_stress_begin =
f_stress.begin(spatial_dimension, spatial_dimension * 2);
std::vector<Real> new_sigmas;
std::vector< Vector<Real> > new_normal_traction;
std::vector<Real> new_delta_c;
// loop over each facet belonging to this material
for (UInt f = 0; f < nb_facet; ++f, ++sigma_lim_it) {
UInt facet = f_filter(f);
// skip facets where check shouldn't be realized
if (!facets_check(facet)) continue;
// compute the effective norm on each quadrature point of the facet
for (UInt q = 0; q < nb_quad_facet; ++q) {
UInt current_quad = facet * nb_quad_facet + q;
const Vector<Real> & normal = normal_begin[current_quad];
const Vector<Real> & tangent = tangent_begin[current_quad];
const Matrix<Real> & facet_stress_it = facet_stress_begin[current_quad];
// compute average stress on the current quadrature point
Matrix<Real> stress_1(facet_stress_it.storage(),
spatial_dimension,
spatial_dimension);
Matrix<Real> stress_2(facet_stress_it.storage() + sp2,
spatial_dimension,
spatial_dimension);
stress_tmp.copy(stress_1);
stress_tmp += stress_2;
stress_tmp /= 2.;
Vector<Real> normal_traction_vec(normal_traction(q));
// compute normal and effective stress
stress_check(q) = computeEffectiveNorm(stress_tmp, normal, tangent,
normal_traction_vec);
}
// verify if the effective stress overcomes the threshold
if (stress_check.mean() > (*sigma_lim_it - tolerance)) {
if (model->isExplicit()){
f_insertion(facet) = true;
// store the new cohesive material parameters for each quadrature point
for (UInt q = 0; q < nb_quad_facet; ++q) {
Real new_sigma = stress_check(q);
Vector<Real> normal_traction_vec(normal_traction(q));
if (spatial_dimension != 3)
normal_traction_vec *= -1.;
new_sigmas.push_back(new_sigma);
new_normal_traction.push_back(normal_traction_vec);
Real new_delta;
// set delta_c in function of G_c or a given delta_c value
if (Math::are_float_equal(delta_c, 0.))
new_delta = 2 * G_c / new_sigma;
else
new_delta = (*sigma_lim_it) / new_sigma * delta_c;
new_delta_c.push_back(new_delta);
}
}else{
Real ratio = stress_check.mean()/(*sigma_lim_it);
if (ratio > max_ratio){
++nn;
max_ratio = ratio;
index_f = f;
index_filter = f_filter(f);
//sigma_c_tmp = *sigma_lim_it;
//sigma_ins_tmp = ins_s_;
}
}
}
}
- StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
- Array<Real> abs_max(comm.getNbProc());
- abs_max(comm.whoAmI()) = max_ratio;
- comm.allGather(abs_max.storage(), 1);
-
- Array<Real>::scalar_iterator it = std::max_element(abs_max.begin(), abs_max.end());
- UInt pos = it - abs_max.begin();
-
- if (pos != comm.whoAmI()) {
- AKANTU_DEBUG_OUT();
- return;
- }
-
/// insertion of only 1 cohesive element in case of implicit approach. The one subjected to the highest stress.
- if (!model->isExplicit() && nn){
- f_insertion(index_filter) = true;
+ if (!model->isExplicit()){
+ StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
+ Array<Real> abs_max(comm.getNbProc());
+ abs_max(comm.whoAmI()) = max_ratio;
+ comm.allGather(abs_max.storage(), 1);
+
+ Array<Real>::scalar_iterator it = std::max_element(abs_max.begin(), abs_max.end());
+ UInt pos = it - abs_max.begin();
+
+ if (pos != comm.whoAmI()) {
+ AKANTU_DEBUG_OUT();
+ return;
+ }
- // Array<Real>::iterator<Matrix<Real> > normal_traction_it =
- // normal_traction.begin_reinterpret(nb_quad_facet, spatial_dimension, nb_facet);
- Array<Real>::const_iterator<Real> sigma_lim_it = sigma_lim.begin();
+ if (nn) {
+ f_insertion(index_filter) = true;
- for (UInt q = 0; q < nb_quad_facet; ++q) {
+ // Array<Real>::iterator<Matrix<Real> > normal_traction_it =
+ // normal_traction.begin_reinterpret(nb_quad_facet, spatial_dimension, nb_facet);
+ Array<Real>::const_iterator<Real> sigma_lim_it = sigma_lim.begin();
- // Vector<Real> ins_s(normal_traction_it[index_f].storage() + q * spatial_dimension,
- // spatial_dimension);
+ for (UInt q = 0; q < nb_quad_facet; ++q) {
- Real new_sigma = (sigma_lim_it[index_f]);
+ // Vector<Real> ins_s(normal_traction_it[index_f].storage() + q * spatial_dimension,
+ // spatial_dimension);
- new_sigmas.push_back(new_sigma);
- new_normal_traction.push_back(0.0);
+ Real new_sigma = (sigma_lim_it[index_f]);
- //// sig_c_eff.push_back(new_sigma);
- // ins_stress.push_back(ins_s);
- // trac_old.push_back(ins_s);
- //// ins_stress.push_back(0.0);
- //// trac_old.push_back(0.0);
+ new_sigmas.push_back(new_sigma);
+ new_normal_traction.push_back(0.0);
- Real new_delta;
+ //// sig_c_eff.push_back(new_sigma);
+ // ins_stress.push_back(ins_s);
+ // trac_old.push_back(ins_s);
+ //// ins_stress.push_back(0.0);
+ //// trac_old.push_back(0.0);
- //set delta_c in function of G_c or a given delta_c value
- if (!Math::are_float_equal(delta_c, 0.))
- new_delta = delta_c;
- else
- new_delta = 2 * G_c / (new_sigma);
+ Real new_delta;
- new_delta_c.push_back(new_delta);
- ////del_c.push_back(new_delta);
+ //set delta_c in function of G_c or a given delta_c value
+ if (!Math::are_float_equal(delta_c, 0.))
+ new_delta = delta_c;
+ else
+ new_delta = 2 * G_c / (new_sigma);
+ new_delta_c.push_back(new_delta);
+ ////del_c.push_back(new_delta);
+
+ }
}
}
-
// update material data for the new elements
UInt old_nb_quad_points = sig_c_eff.getSize();
UInt new_nb_quad_points = new_sigmas.size();
sig_c_eff.resize(old_nb_quad_points + new_nb_quad_points);
ins_stress.resize(old_nb_quad_points + new_nb_quad_points);
trac_old.resize(old_nb_quad_points + new_nb_quad_points);
del_c.resize(old_nb_quad_points + new_nb_quad_points);
for (UInt q = 0; q < new_nb_quad_points; ++q) {
sig_c_eff(old_nb_quad_points + q) = new_sigmas[q];
del_c(old_nb_quad_points + q) = new_delta_c[q];
for (UInt dim = 0; dim < spatial_dimension; ++dim) {
ins_stress(old_nb_quad_points + q, dim) = new_normal_traction[q](dim);
trac_old(old_nb_quad_points + q, dim) = new_normal_traction[q](dim);
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
inline Real MaterialCohesiveLinear<spatial_dimension>::computeEffectiveNorm(const Matrix<Real> & stress,
const Vector<Real> & normal,
const Vector<Real> & tangent,
Vector<Real> & normal_traction) {
normal_traction.mul<false>(stress, normal);
Real normal_contrib = normal_traction.dot(normal);
/// in 3D tangential components must be summed
Real tangent_contrib = 0;
if (spatial_dimension == 2) {
Real tangent_contrib_tmp = normal_traction.dot(tangent);
tangent_contrib += tangent_contrib_tmp * tangent_contrib_tmp;
}
else if (spatial_dimension == 3) {
for (UInt s = 0; s < spatial_dimension - 1; ++s) {
const Vector<Real> tangent_v(tangent.storage() + s * spatial_dimension,
spatial_dimension);
Real tangent_contrib_tmp = normal_traction.dot(tangent_v);
tangent_contrib += tangent_contrib_tmp * tangent_contrib_tmp;
}
}
tangent_contrib = std::sqrt(tangent_contrib);
normal_contrib = std::max(0., normal_contrib);
return std::sqrt(normal_contrib * normal_contrib + tangent_contrib * tangent_contrib * beta2_inv);
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialCohesiveLinear<spatial_dimension>::computeTraction(const Array<Real> & normal,
ElementType el_type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
/// define iterators
Array<Real>::iterator< Vector<Real> > traction_it =
tractions(el_type, ghost_type).begin(spatial_dimension);
Array<Real>::iterator< Vector<Real> > opening_it =
opening(el_type, ghost_type).begin(spatial_dimension);
Array<Real>::iterator< Vector<Real> > contact_traction_it =
contact_tractions(el_type, ghost_type).begin(spatial_dimension);
Array<Real>::iterator< Vector<Real> > contact_opening_it =
contact_opening(el_type, ghost_type).begin(spatial_dimension);
Array<Real>::const_iterator< Vector<Real> > normal_it =
normal.begin(spatial_dimension);
Array<Real>::iterator< Vector<Real> >traction_end =
tractions(el_type, ghost_type).end(spatial_dimension);
Array<Real>::iterator<Real>sigma_c_it =
sigma_c_eff(el_type, ghost_type).begin();
Array<Real>::iterator<Real>delta_max_it =
delta_max(el_type, ghost_type).begin();
Array<Real>::iterator<Real>delta_max_prev_it =
delta_max.previous(el_type, ghost_type).begin();
Array<Real>::iterator<Real>delta_c_it =
delta_c_eff(el_type, ghost_type).begin();
Array<Real>::iterator<Real>damage_it =
damage(el_type, ghost_type).begin();
Array<Real>::iterator<Vector<Real> > insertion_stress_it =
insertion_stress(el_type, ghost_type).begin(spatial_dimension);
Real * memory_space = new Real[2*spatial_dimension];
Vector<Real> normal_opening(memory_space, spatial_dimension);
Vector<Real> tangential_opening(memory_space + spatial_dimension,
spatial_dimension);
/// loop on each quadrature point
for (; traction_it != traction_end;
++traction_it, ++opening_it, ++normal_it, ++sigma_c_it,
++delta_max_it, ++delta_c_it, ++damage_it, ++contact_traction_it,
++insertion_stress_it, ++contact_opening_it, ++delta_max_prev_it) {
if (!model->isExplicit())
*delta_max_it = *delta_max_prev_it;
/// compute normal and tangential opening vectors
Real normal_opening_norm = opening_it->dot(*normal_it);
normal_opening = (*normal_it);
normal_opening *= normal_opening_norm;
tangential_opening = *opening_it;
tangential_opening -= normal_opening;
Real tangential_opening_norm = tangential_opening.norm();
/**
* compute effective opening displacement
* @f$ \delta = \sqrt{
* \frac{\beta^2}{\kappa^2} \Delta_t^2 + \Delta_n^2 } @f$
*/
Real delta = tangential_opening_norm * tangential_opening_norm * beta2_kappa2;
bool penetration = normal_opening_norm < -Math::getTolerance();
if (penetration) {
/// use penalty coefficient in case of penetration
*contact_traction_it = normal_opening;
*contact_traction_it *= penalty;
*contact_opening_it = normal_opening;
/// don't consider penetration contribution for delta
*opening_it = tangential_opening;
normal_opening.clear();
}
else {
delta += normal_opening_norm * normal_opening_norm;
contact_traction_it->clear();
contact_opening_it->clear();
}
delta = std::sqrt(delta);
/// update maximum displacement and damage
*delta_max_it = std::max(*delta_max_it, delta);
*damage_it = std::min(*delta_max_it / *delta_c_it, 1.);
/**
* Compute traction @f$ \mathbf{T} = \left(
* \frac{\beta^2}{\kappa} \Delta_t \mathbf{t} + \Delta_n
* \mathbf{n} \right) \frac{\sigma_c}{\delta} \left( 1-
* \frac{\delta}{\delta_c} \right)@f$
*/
if (Math::are_float_equal(*damage_it, 1.))
traction_it->clear();
else if (Math::are_float_equal(*damage_it, 0.)) {
if (penetration)
traction_it->clear();
else
*traction_it = *insertion_stress_it;
}
else {
*traction_it = tangential_opening;
*traction_it *= beta2_kappa;
*traction_it += normal_opening;
AKANTU_DEBUG_ASSERT(*delta_max_it != 0.,
"Division by zero, tolerance might be too low");
*traction_it *= *sigma_c_it / *delta_max_it * (1. - *damage_it);
}
}
delete [] memory_space;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialCohesiveLinear<spatial_dimension>::computeTangentTraction(const ElementType & el_type,
Array<Real> & tangent_matrix,
const Array<Real> & normal,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
/// define iterators
Array<Real>::matrix_iterator tangent_it =
tangent_matrix.begin(spatial_dimension, spatial_dimension);
Array<Real>::matrix_iterator tangent_end =
tangent_matrix.end(spatial_dimension, spatial_dimension);
Array<Real>::const_vector_iterator normal_it =
normal.begin(spatial_dimension);
Array<Real>::vector_iterator opening_it =
opening(el_type, ghost_type).begin(spatial_dimension);
Array<Real>::iterator<Real>delta_max_it =
delta_max.previous(el_type, ghost_type).begin();
Array<Real>::iterator<Real>sigma_c_it =
sigma_c_eff(el_type, ghost_type).begin();
Array<Real>::iterator<Real>delta_c_it =
delta_c_eff(el_type, ghost_type).begin();
Array<Real>::iterator<Real>damage_it =
damage(el_type, ghost_type).begin();
Array<Real>::iterator< Vector<Real> > contact_opening_it =
contact_opening(el_type, ghost_type).begin(spatial_dimension);
Vector<Real> normal_opening(spatial_dimension);
Vector<Real> tangential_opening(spatial_dimension);
for (; tangent_it != tangent_end;
++tangent_it, ++normal_it, ++opening_it, ++ delta_max_it,
++sigma_c_it, ++delta_c_it, ++damage_it, ++contact_opening_it) {
/// compute normal and tangential opening vectors
*opening_it += *contact_opening_it;
Real normal_opening_norm = opening_it->dot(*normal_it);
normal_opening = (*normal_it);
normal_opening *= normal_opening_norm;
tangential_opening = *opening_it;
tangential_opening -= normal_opening;
Real tangential_opening_norm = tangential_opening.norm();
bool penetration = normal_opening_norm < -Math::getTolerance();
Real derivative = 0;
Real t = 0;
Real delta = tangential_opening_norm * tangential_opening_norm * beta2_kappa2;
delta += normal_opening_norm * normal_opening_norm;
delta = std::sqrt(delta);
if (delta < Math::getTolerance())
delta = (*delta_c_it)/1000.; //0.0000001;
if (normal_opening_norm >= 0.0){
if (delta >= *delta_max_it){
derivative = -*sigma_c_it/(delta * delta);
t = *sigma_c_it * (1 - delta / *delta_c_it);
} else if (delta < *delta_max_it){
Real tmax = *sigma_c_it * (1 - *delta_max_it / *delta_c_it);
t = tmax / *delta_max_it * delta;
}
}
Matrix<Real> n_outer_n(spatial_dimension, spatial_dimension);
n_outer_n.outerProduct(*normal_it, *normal_it);
if (penetration){
///don't consider penetration contribution for delta
*opening_it = tangential_opening;
*tangent_it += n_outer_n;
*tangent_it *= penalty;
}
Matrix<Real> I(spatial_dimension, spatial_dimension);
I.eye(beta2_kappa);
Matrix<Real> nn(n_outer_n);
nn *= (1 - beta2_kappa);
nn += I;
nn *= t/delta;
Vector<Real> t_tilde(normal_opening);
t_tilde *= (1 - beta2_kappa2);
Vector<Real> mm(*opening_it);
mm *= beta2_kappa2;
t_tilde += mm;
Vector<Real> t_hat(normal_opening);
t_hat += beta2_kappa * tangential_opening;
Matrix<Real> prov(spatial_dimension, spatial_dimension);
prov.outerProduct(t_hat, t_tilde);
prov *= derivative/delta;
prov += nn;
*tangent_it += prov;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
INSTANSIATE_MATERIAL(MaterialCohesiveLinear);
__END_AKANTU__
diff --git a/src/model/solid_mechanics/materials/material_embedded/embedded_internal_field.hh b/src/model/solid_mechanics/materials/material_embedded/embedded_internal_field.hh
index dd1d2a83f..2692bcc3a 100644
--- a/src/model/solid_mechanics/materials/material_embedded/embedded_internal_field.hh
+++ b/src/model/solid_mechanics/materials/material_embedded/embedded_internal_field.hh
@@ -1,74 +1,81 @@
/**
* @file embedded_internal_field.hh
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @date creation: Thu Mar 19 2015
* @date last modification: Thu Mar 19 2015
*
* @brief Embedded Material internal properties
*
* @section LICENSE
*
* Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "internal_field.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_EMBEDDED_INTERNAL_FIELD_HH__
#define __AKANTU_EMBEDDED_INTERNAL_FIELD_HH__
__BEGIN_AKANTU__
class Material;
class FEEngine;
template<typename T>
class EmbeddedInternalField : public InternalField<T> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
EmbeddedInternalField(const ID & id, Material & material):
InternalField<T>(id,
material,
material.getModel().getFEEngine("EmbeddedInterfaceFEEngine"),
material.getElementFilter())
{
this->spatial_dimension = 1;
}
EmbeddedInternalField(const ID & id, const EmbeddedInternalField & other):
InternalField<T>(id, other)
{
this->spatial_dimension = 1;
}
void operator=(const EmbeddedInternalField & other) {
InternalField<T>::operator=(other);
this->spatial_dimension = 1;
}
};
+template<>
+inline void ParsableParamTyped< EmbeddedInternalField<Real> >::parseParam(const ParserParameter & in_param) {
+ ParsableParam::parseParam(in_param);
+ Real r = in_param;
+ param.setDefaultValue(r);
+}
+
__END_AKANTU__
#endif // __AKANTU_EMBEDDED_INTERNAL_FIELD_HH__
diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc
index 347ee03f2..2b160ef21 100644
--- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc
+++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc
@@ -1,673 +1,676 @@
/**
* @file material_reinforcement.cc
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @date creation: Thu Mar 12 2015
* @date last modification: Thu Mar 12 2015
*
* @brief Reinforcement material
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_voigthelper.hh"
#include "material_reinforcement.hh"
__BEGIN_AKANTU__
template<UInt dim>
MaterialReinforcement<dim>::MaterialReinforcement(SolidMechanicsModel & model, const ID & id):
Material(model, id),
model(NULL),
gradu("gradu_embedded", *this),
stress("stress_embedded", *this),
directing_cosines("directing_cosines", *this),
+ pre_stress("pre_stress", *this),
area(1.0),
shape_derivatives()
{
this->model = dynamic_cast<EmbeddedInterfaceModel *>(&model);
AKANTU_DEBUG_ASSERT(this->model != NULL, "MaterialReinforcement needs an EmbeddedInterfaceModel");
this->model->getInterfaceMesh().initElementTypeMapArray(element_filter, 1, 1,
false, _ek_regular);
this->registerParam("area", area, _pat_parsable | _pat_modifiable, "Reinforcement cross-sectional area");
+ this->registerParam("pre_stress", pre_stress, _pat_parsable | _pat_modifiable, "Uniform pre-stress");
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
MaterialReinforcement<dim>::~MaterialReinforcement() {
AKANTU_DEBUG_IN();
ElementTypeMap<ElementTypeMapArray<Real> *>::type_iterator it = shape_derivatives.firstType();
ElementTypeMap<ElementTypeMapArray<Real> *>::type_iterator end = shape_derivatives.lastType();
for (; it != end ; ++it) {
delete shape_derivatives(*it, _not_ghost);
delete shape_derivatives(*it, _ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void MaterialReinforcement<dim>::initMaterial() {
Material::initMaterial();
gradu.initialize(dim * dim);
stress.initialize(dim * dim);
+ pre_stress.initialize(1);
/// We initialise the stuff that is not going to change during the simulation
this->allocBackgroundShapeDerivatives();
this->initBackgroundShapeDerivatives();
this->initDirectingCosines();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void MaterialReinforcement<dim>::allocBackgroundShapeDerivatives() {
AKANTU_DEBUG_IN();
Mesh & interface_mesh = model->getInterfaceMesh();
Mesh & mesh = model->getMesh();
ghost_type_t::iterator int_ghost_it = ghost_type_t::begin();
// Loop over interface ghosts
for (; int_ghost_it != ghost_type_t::end() ; ++int_ghost_it) {
Mesh::type_iterator interface_type_it = interface_mesh.firstType();
Mesh::type_iterator interface_type_end = interface_mesh.lastType();
for (; interface_type_it != interface_type_end ; ++interface_type_it) {
Mesh::type_iterator background_type_it = mesh.firstType(dim, *int_ghost_it);
Mesh::type_iterator background_type_end = mesh.lastType(dim, *int_ghost_it);
for (; background_type_it != background_type_end ; ++background_type_it) {
const ElementType & int_type = *interface_type_it;
const ElementType & back_type = *background_type_it;
const GhostType & int_ghost = *int_ghost_it;
std::string shaped_id = "embedded_shape_derivatives";
if (int_ghost == _ghost) shaped_id += ":ghost";
ElementTypeMapArray<Real> * shaped_etma = new ElementTypeMapArray<Real>(shaped_id, this->name);
UInt nb_points = Mesh::getNbNodesPerElement(back_type);
UInt nb_quad_points = model->getFEEngine("EmbeddedInterfaceFEEngine").getNbQuadraturePoints(int_type);
UInt nb_elements = element_filter(int_type, int_ghost).getSize();
shaped_etma->alloc(nb_elements * nb_quad_points,
dim * nb_points,
back_type);
shape_derivatives(shaped_etma, int_type, int_ghost);
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void MaterialReinforcement<dim>::initBackgroundShapeDerivatives() {
AKANTU_DEBUG_IN();
Mesh & mesh = model->getMesh();
Mesh::type_iterator type_it = mesh.firstType(dim, _not_ghost);
Mesh::type_iterator type_end = mesh.lastType(dim, _not_ghost);
for (; type_it != type_end ; ++type_it) {
computeBackgroundShapeDerivatives(*type_it, _not_ghost);
//computeBackgroundShapeDerivatives(*type_it, _ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void MaterialReinforcement<dim>::initDirectingCosines() {
AKANTU_DEBUG_IN();
Mesh & mesh = model->getInterfaceMesh();
Mesh::type_iterator type_it = mesh.firstType(1, _not_ghost);
Mesh::type_iterator type_end = mesh.lastType(1, _not_ghost);
const UInt voigt_size = getTangentStiffnessVoigtSize(dim);
directing_cosines.initialize(voigt_size * voigt_size);
for (; type_it != type_end ; ++type_it) {
computeDirectingCosines(*type_it, _not_ghost);
computeDirectingCosines(*type_it, _ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void MaterialReinforcement<dim>::assembleStiffnessMatrix(GhostType ghost_type) {
AKANTU_DEBUG_IN();
Mesh & interface_mesh = model->getInterfaceMesh();
Mesh::type_iterator type_it = interface_mesh.firstType();
Mesh::type_iterator type_end = interface_mesh.lastType();
for (; type_it != type_end ; ++type_it) {
assembleStiffnessMatrix(*type_it, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void MaterialReinforcement<dim>::updateResidual(GhostType ghost_type) {
AKANTU_DEBUG_IN();
computeAllStresses(ghost_type);
assembleResidual(ghost_type);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void MaterialReinforcement<dim>::assembleResidual(GhostType ghost_type) {
AKANTU_DEBUG_IN();
Mesh & interface_mesh = model->getInterfaceMesh();
Mesh::type_iterator type_it = interface_mesh.firstType();
Mesh::type_iterator type_end = interface_mesh.lastType();
for (; type_it != type_end ; ++type_it) {
assembleResidual(*type_it, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void MaterialReinforcement<dim>::computeGradU(const ElementType & type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
Array<UInt> & elem_filter = element_filter(type, ghost_type);
UInt nb_element = elem_filter.getSize();
UInt nb_quad_points = model->getFEEngine("EmbeddedInterfaceFEEngine").getNbQuadraturePoints(type);
Array<Real> & gradu_vec = gradu(type, ghost_type);
Mesh::type_iterator back_it = model->getMesh().firstType(dim, ghost_type);
Mesh::type_iterator back_end = model->getMesh().lastType(dim, ghost_type);
for (; back_it != back_end ; ++back_it) {
UInt nodes_per_background_e = Mesh::getNbNodesPerElement(*back_it);
Array<Real> * shapesd_filtered = new Array<Real>(nb_element, dim * nodes_per_background_e, "shapesd_filtered");
FEEngine::filterElementalData(model->getInterfaceMesh(),
shape_derivatives(type, ghost_type)->operator()(*back_it, ghost_type),
*shapesd_filtered,
type, ghost_type, elem_filter);
Array<UInt> * background_filter = new Array<UInt>(nb_element, 1, "background_filter");
filterInterfaceBackgroundElements(*background_filter, *back_it, type, ghost_type, ghost_type);
Array<Real> * disp_per_element = new Array<Real>(0, dim * nodes_per_background_e, "disp_elem");
FEEngine::extractNodalToElementField(model->getMesh(),
model->getDisplacement(),
*disp_per_element,
*back_it, ghost_type, *background_filter);
Array<Real>::matrix_iterator disp_it = disp_per_element->begin(dim, nodes_per_background_e);
Array<Real>::matrix_iterator disp_end = disp_per_element->end(dim, nodes_per_background_e);
Array<Real>::matrix_iterator shapes_it = shapesd_filtered->begin(dim, nodes_per_background_e);
Array<Real>::matrix_iterator grad_u_it = gradu_vec.begin(dim, dim);
for (; disp_it != disp_end ; ++disp_it) {
for (UInt i = 0; i < nb_quad_points; i++, ++shapes_it, ++grad_u_it) {
Matrix<Real> & B = *shapes_it;
Matrix<Real> & du = *grad_u_it;
Matrix<Real> & u = *disp_it;
du.mul<false, true>(u, B);
}
}
delete shapesd_filtered;
delete background_filter;
delete disp_per_element;
}
AKANTU_DEBUG_OUT();
}
template<UInt dim>
void MaterialReinforcement<dim>::computeAllStresses(GhostType ghost_type) {
AKANTU_DEBUG_IN();
Mesh::type_iterator it = model->getInterfaceMesh().firstType();
Mesh::type_iterator last_type = model->getInterfaceMesh().lastType();
for(; it != last_type; ++it) {
computeGradU(*it, ghost_type);
computeStress(*it, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void MaterialReinforcement<dim>::assembleResidual(const ElementType & type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
Mesh & mesh = model->getMesh();
Mesh::type_iterator type_it = mesh.firstType(dim, ghost_type);
Mesh::type_iterator type_end = mesh.lastType(dim, ghost_type);
for (; type_it != type_end ; ++type_it) {
assembleResidual(type, *type_it, ghost_type, _not_ghost);
//assembleResidual(type, *type_it, ghost_type, _ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void MaterialReinforcement<dim>::assembleResidual(const ElementType & interface_type,
const ElementType & background_type,
GhostType interface_ghost,
GhostType background_ghost) {
AKANTU_DEBUG_IN();
UInt voigt_size = getTangentStiffnessVoigtSize(dim);
Array<Real> & residual = const_cast<Array<Real> &>(model->getResidual());
FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine");
FEEngine & background_engine = model->getFEEngine();
Array<UInt> & elem_filter = element_filter(interface_type, interface_ghost);
UInt nodes_per_background_e = Mesh::getNbNodesPerElement(background_type);
UInt nb_quadrature_points = interface_engine.getNbQuadraturePoints(interface_type, interface_ghost);
UInt nb_element = elem_filter.getSize();
UInt back_dof = dim * nodes_per_background_e;
Array<Real> & shapesd = shape_derivatives(interface_type, interface_ghost)->operator()(background_type, background_ghost);
Array<Real> * shapesd_filtered = new Array<Real>(nb_element * nb_quadrature_points,
back_dof,
"background_shapesd");
FEEngine::filterElementalData(model->getInterfaceMesh(), shapesd, *shapesd_filtered,
interface_type, interface_ghost, elem_filter);
Array<Real> * integrant = new Array<Real>(nb_quadrature_points * nb_element,
back_dof,
"integrant");
Array<Real>::vector_iterator integrant_it =
integrant->begin(back_dof);
Array<Real>::vector_iterator integrant_end =
integrant->end(back_dof);
Array<Real>::matrix_iterator B_it =
shapesd_filtered->begin(dim, nodes_per_background_e);
Array<Real>::matrix_iterator C_it =
directing_cosines(interface_type, interface_ghost).begin(voigt_size, voigt_size);
Array<Real>::matrix_iterator sigma_it =
stress(interface_type, interface_ghost).begin(dim, dim);
Vector<Real> sigma(voigt_size);
Matrix<Real> Bvoigt(voigt_size, back_dof);
Vector<Real> Ct_sigma(voigt_size);
for (; integrant_it != integrant_end ; ++integrant_it,
++B_it,
++C_it,
++sigma_it) {
VoigtHelper<dim>::transferBMatrixToSymVoigtBMatrix(*B_it, Bvoigt, nodes_per_background_e);
Matrix<Real> & C = *C_it;
Vector<Real> & BtCt_sigma = *integrant_it;
stressTensorToVoigtVector(*sigma_it, sigma);
Ct_sigma.mul<true>(C, sigma);
BtCt_sigma.mul<true>(Bvoigt, Ct_sigma);
BtCt_sigma *= area;
}
delete shapesd_filtered;
Array<Real> * residual_interface = new Array<Real>(nb_element, back_dof, "residual_interface");
interface_engine.integrate(*integrant,
*residual_interface,
back_dof,
interface_type, interface_ghost,
elem_filter);
delete integrant;
Array<UInt> * background_filter = new Array<UInt>(nb_element, 1, "background_filter");
filterInterfaceBackgroundElements(*background_filter,
background_type, interface_type,
background_ghost, interface_ghost);
background_engine.assembleArray(*residual_interface, residual,
model->getDOFSynchronizer().getLocalDOFEquationNumbers(),
dim, background_type, background_ghost, *background_filter, -1.0);
delete residual_interface;
delete background_filter;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void MaterialReinforcement<dim>::filterInterfaceBackgroundElements(Array<UInt> & filter,
const ElementType & type,
const ElementType & interface_type,
GhostType ghost_type,
GhostType interface_ghost_type) {
AKANTU_DEBUG_IN();
filter.resize(0);
filter.clear();
Array<Element> & elements = model->getInterfaceAssociatedElements(interface_type, interface_ghost_type);
Array<UInt> & elem_filter = element_filter(interface_type, interface_ghost_type);
Array<UInt>::scalar_iterator
filter_it = elem_filter.begin(),
filter_end = elem_filter.end();
for (; filter_it != filter_end ; ++filter_it) {
Element & elem = elements(*filter_it);
if (elem.type == type) filter.push_back(elem.element);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void MaterialReinforcement<dim>::computeDirectingCosines(const ElementType & type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
Mesh & interface_mesh = this->model->getInterfaceMesh();
const UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
const UInt steel_dof = dim * nb_nodes_per_element;
const UInt voigt_size = getTangentStiffnessVoigtSize(dim);
const UInt nb_quad_points =
model->getFEEngine("EmbeddedInterfaceFEEngine").getNbQuadraturePoints(type, ghost_type);
Array<Real> node_coordinates(this->element_filter(type, ghost_type).getSize(), steel_dof);
this->model->getFEEngine().template extractNodalToElementField<Real>(interface_mesh,
interface_mesh.getNodes(),
node_coordinates,
type,
ghost_type,
this->element_filter(type, ghost_type));
Array<Real>::matrix_iterator
directing_cosines_it = directing_cosines(type, ghost_type).begin(voigt_size, voigt_size);
Array<Real>::matrix_iterator node_coordinates_it = node_coordinates.begin(dim, nb_nodes_per_element);
Array<Real>::matrix_iterator node_coordinates_end = node_coordinates.end(dim, nb_nodes_per_element);
for (; node_coordinates_it != node_coordinates_end ; ++node_coordinates_it) {
for (UInt i = 0 ; i < nb_quad_points ; i++, ++directing_cosines_it) {
Matrix<Real> & nodes = *node_coordinates_it;
Matrix<Real> & cosines = *directing_cosines_it;
computeDirectingCosinesOnQuad(nodes, cosines);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void MaterialReinforcement<dim>::assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
Mesh & mesh = model->getMesh();
Mesh::type_iterator type_it = mesh.firstType(dim, ghost_type);
Mesh::type_iterator type_end = mesh.lastType(dim, ghost_type);
for (; type_it != type_end ; ++type_it) {
assembleStiffnessMatrix(type, *type_it, ghost_type, _not_ghost);
//assembleStiffnessMatrix(type, *type_it, ghost_type, _ghost);
}
AKANTU_DEBUG_OUT();
}
template<UInt dim>
void MaterialReinforcement<dim>::assembleStiffnessMatrix(const ElementType & interface_type,
const ElementType & background_type,
GhostType interface_ghost,
GhostType background_ghost) {
AKANTU_DEBUG_IN();
UInt voigt_size = getTangentStiffnessVoigtSize(dim);
SparseMatrix & K = const_cast<SparseMatrix &>(model->getStiffnessMatrix());
FEEngine & background_engine = model->getFEEngine();
FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine");
Array<UInt> & elem_filter = element_filter(interface_type, interface_ghost);
Array<Real> & grad_u = gradu(interface_type, interface_ghost);
UInt nb_element = elem_filter.getSize();
UInt nodes_per_background_e = Mesh::getNbNodesPerElement(background_type);
UInt nb_quadrature_points = interface_engine.getNbQuadraturePoints(interface_type, interface_ghost);
UInt back_dof = dim * nodes_per_background_e;
UInt integrant_size = back_dof;
grad_u.resize(nb_quadrature_points * nb_element);
//need function model->getInterfaceDisplacement()
/*interface_engine.gradientOnQuadraturePoints(model->getInterfaceDisplacement(),
grad_u, dim, interface_type, interface_ghost,
elem_filter);*/
Array<Real> * tangent_moduli = new Array<Real>(nb_element * nb_quadrature_points,
1, "interface_tangent_moduli");
tangent_moduli->clear();
computeTangentModuli(interface_type, *tangent_moduli, interface_ghost);
Array<Real> & shapesd = shape_derivatives(interface_type, interface_ghost)->operator()(background_type, background_ghost);
Array<Real> * shapesd_filtered = new Array<Real>(nb_element * nb_quadrature_points,
back_dof,
"background_shapesd");
FEEngine::filterElementalData(model->getInterfaceMesh(), shapesd, *shapesd_filtered,
interface_type, interface_ghost, elem_filter);
Array<Real> * integrant = new Array<Real>(nb_element * nb_quadrature_points,
integrant_size * integrant_size,
"B^t*C^t*D*C*B");
integrant->clear();
/// Temporary matrices for integrant product
Matrix<Real> Bvoigt(voigt_size, back_dof);
Matrix<Real> DC(voigt_size, voigt_size);
Matrix<Real> DCB(voigt_size, back_dof);
Matrix<Real> CtDCB(voigt_size, back_dof);
Array<Real>::scalar_iterator D_it = tangent_moduli->begin();
Array<Real>::scalar_iterator D_end = tangent_moduli->end();
Array<Real>::matrix_iterator C_it =
directing_cosines(interface_type, interface_ghost).begin(voigt_size, voigt_size);
Array<Real>::matrix_iterator B_it =
shapesd_filtered->begin(dim, nodes_per_background_e);
Array<Real>::matrix_iterator integrant_it =
integrant->begin(integrant_size, integrant_size);
for (; D_it != D_end ; ++D_it, ++C_it, ++B_it, ++integrant_it) {
Real & D = *D_it;
Matrix<Real> & C = *C_it;
Matrix<Real> & B = *B_it;
Matrix<Real> & BtCtDCB = *integrant_it;
VoigtHelper<dim>::transferBMatrixToSymVoigtBMatrix(B, Bvoigt, nodes_per_background_e);
DC.clear();
DC(0, 0) = D * area;
DC *= C;
DCB.mul<false, false>(DC, Bvoigt);
CtDCB.mul<true, false>(C, DCB);
BtCtDCB.mul<true, false>(Bvoigt, CtDCB);
}
delete tangent_moduli;
delete shapesd_filtered;
Array<Real> * K_interface = new Array<Real>(nb_element, integrant_size * integrant_size, "K_interface");
interface_engine.integrate(*integrant, *K_interface,
integrant_size * integrant_size,
interface_type, interface_ghost,
elem_filter);
delete integrant;
Array<UInt> * background_filter = new Array<UInt>(nb_element, 1, "background_filter");
filterInterfaceBackgroundElements(*background_filter,
background_type, interface_type,
background_ghost, interface_ghost);
background_engine.assembleMatrix(*K_interface, K, dim, background_type, background_ghost, *background_filter);
delete K_interface;
delete background_filter;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/// In this function, type and ghost type refer to background elements
template<UInt dim>
void MaterialReinforcement<dim>::computeBackgroundShapeDerivatives(const ElementType & type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
Mesh & interface_mesh = model->getInterfaceMesh();
FEEngine & engine = model->getFEEngine();
FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine");
Mesh::type_iterator interface_type = interface_mesh.firstType();
Mesh::type_iterator interface_last = interface_mesh.lastType();
for (; interface_type != interface_last ; ++interface_type) {
Array<UInt> & filter = element_filter(*interface_type, ghost_type);
const UInt nb_elements = filter.getSize();
const UInt nb_nodes = Mesh::getNbNodesPerElement(type);
const UInt nb_quad_per_element = interface_engine.getNbQuadraturePoints(*interface_type);
Array<Real> quad_pos(nb_quad_per_element * nb_elements, dim, "interface_quad_points");
quad_pos.resize(nb_quad_per_element * nb_elements);
interface_engine.interpolateOnQuadraturePoints(interface_mesh.getNodes(),
quad_pos, dim, *interface_type,
ghost_type, filter);
Array<Real> & background_shapesd = shape_derivatives(*interface_type, ghost_type)->operator()(type, ghost_type);
background_shapesd.clear();
Array<UInt> * background_elements = new Array<UInt>(nb_elements, 1, "computeBackgroundShapeDerivatives:background_filter");
filterInterfaceBackgroundElements(*background_elements, type, *interface_type, ghost_type, ghost_type);
Array<UInt>::scalar_iterator
back_it = background_elements->begin(),
back_end = background_elements->end();
Array<Real>::matrix_iterator shapesd_it = background_shapesd.begin(dim, nb_nodes);
Array<Real>::vector_iterator quad_pos_it = quad_pos.begin(dim);
for (; back_it != back_end ; ++back_it) {
for (UInt i = 0 ; i < nb_quad_per_element ; i++, ++shapesd_it, ++quad_pos_it)
engine.computeShapeDerivatives(*quad_pos_it, *back_it, type, *shapesd_it, ghost_type);
}
delete background_elements;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
INSTANSIATE_MATERIAL(MaterialReinforcement);
/* -------------------------------------------------------------------------- */
__END_AKANTU__
diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh
index 0ccfc3409..1283badbe 100644
--- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh
+++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh
@@ -1,166 +1,169 @@
/**
* @file material_reinforcement.hh
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @date creation: Thu Mar 12 2015
* @date last modification: Thu Mar 12 2015
*
* @brief Reinforcement material
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_REINFORCEMENT_HH__
#define __AKANTU_MATERIAL_REINFORCEMENT_HH__
#include "aka_common.hh"
#include "material.hh"
#include "embedded_interface_model.hh"
#include "embedded_internal_field.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
template<UInt dim>
class MaterialReinforcement : virtual public Material {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
/// Constructor
MaterialReinforcement(SolidMechanicsModel & model, const ID & id = "");
/// Destructor
virtual ~MaterialReinforcement();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// Init the material
virtual void initMaterial();
/// Init the background shape derivatives
void initBackgroundShapeDerivatives();
/// Init the cosine matrices
void initDirectingCosines();
/// Assemble stiffness matrix
virtual void assembleStiffnessMatrix(GhostType ghost_type);
/// Update the residual
virtual void updateResidual(GhostType ghost_type = _not_ghost);
/// Assembled the residual
virtual void assembleResidual(GhostType ghost_type);
/// Compute all the stresses !
virtual void computeAllStresses(GhostType ghost_type);
/// Compute the stiffness parameter for elements of a type
virtual void computeTangentModuli(const ElementType & type,
Array<Real> & tangent,
GhostType ghost_type) = 0;
/* ------------------------------------------------------------------------ */
/* Protected methods */
/* ------------------------------------------------------------------------ */
protected:
/// Allocate the background shape derivatives
void allocBackgroundShapeDerivatives();
/// Compute the directing cosines matrix for one element type
void computeDirectingCosines(const ElementType & type, GhostType ghost_type);
/// Compute the directing cosines matrix on quadrature points
inline void computeDirectingCosinesOnQuad(const Matrix<Real> & nodes,
Matrix<Real> & cosines);
/// Assemble the stiffness matrix for an element type (typically _segment_2)
void assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type);
/// Assemble the stiffness matrix for background & interface types
void assembleStiffnessMatrix(const ElementType & interface_type,
const ElementType & background_type,
GhostType interface_ghost,
GhostType background_ghost);
/// Compute the background shape derivatives for a type (typically _triangle_3 / _tetrahedron_4)
void computeBackgroundShapeDerivatives(const ElementType & type, GhostType ghost_type);
/// Filter elements crossed by interface of a type
void filterInterfaceBackgroundElements(Array<UInt> & filter,
const ElementType & type,
const ElementType & interface_type,
GhostType ghost_type,
GhostType interface_ghost_type);
/// Assemble the residual of one type of element (typically _segment_2)
void assembleResidual(const ElementType & type, GhostType ghost_type);
/// Assemble the residual for a pair of elements
void assembleResidual(const ElementType & interface_type,
const ElementType & background_type,
GhostType interface_ghost,
GhostType background_ghost);
/// TODO figure out why voigt size is 4 in 2D
inline void stressTensorToVoigtVector(const Matrix<Real> & tensor, Vector<Real> & vector);
inline void strainTensorToVoigtVector(const Matrix<Real> & tensor, Vector<Real> & vector);
/// Compute gradu on the interface quadrature points
virtual void computeGradU(const ElementType & type, GhostType ghost_type);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// Embedded model
EmbeddedInterfaceModel * model;
/// grad_u
EmbeddedInternalField<Real> gradu;
/// stress
EmbeddedInternalField<Real> stress;
/// C matrix on quad
EmbeddedInternalField<Real> directing_cosines;
+ /// Prestress on quad
+ EmbeddedInternalField<Real> pre_stress;
+
/// Cross-sectional area
Real area;
/// Background mesh shape derivatives
ElementTypeMap< ElementTypeMapArray<Real> * > shape_derivatives;
};
#include "material_reinforcement_inline_impl.cc"
__END_AKANTU__
#endif // __AKANTU_MATERIAL_REINFORCEMENT_HH__
diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_template_inline_impl.cc b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_template_inline_impl.cc
index 5fa2de99b..98212fa47 100644
--- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_template_inline_impl.cc
+++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_template_inline_impl.cc
@@ -1,178 +1,180 @@
/**
* @file material_reinforcement_template.cc
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @date creation: Mon Mar 16 2015
* @date last modification: Mon Mar 16 2015
*
* @brief Reinforcement material templated with constitutive law
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
template<UInt dim, class ConstLaw>
MaterialReinforcementTemplate<dim, ConstLaw>::MaterialReinforcementTemplate(SolidMechanicsModel & model,
const ID & id):
Material(model, id),
MaterialReinforcement<dim>(model, id),
ConstLaw(model, id)
{}
/* -------------------------------------------------------------------------- */
template<UInt dim, class ConstLaw>
MaterialReinforcementTemplate<dim, ConstLaw>::~MaterialReinforcementTemplate()
{}
/* -------------------------------------------------------------------------- */
template<UInt dim, class ConstLaw>
void MaterialReinforcementTemplate<dim, ConstLaw>::initMaterial() {
MaterialReinforcement<dim>::initMaterial();
// Gotta change the dimension from here onwards
this->ConstLaw::spatial_dimension = 1;
ConstLaw::initMaterial();
this->ConstLaw::gradu.free();
this->ConstLaw::stress.free();
this->ConstLaw::delta_T.free();
this->ConstLaw::sigma_th.free();
this->ConstLaw::potential_energy.free();
Mesh::type_iterator type_it = this->MaterialReinforcement<dim>::model->getInterfaceMesh().firstType();
Mesh::type_iterator type_end = this->MaterialReinforcement<dim>::model->getInterfaceMesh().lastType();
// Reshape the ConstLaw internal fields
for (; type_it != type_end ; ++type_it) {
UInt nb_elements = this->MaterialReinforcement<dim>::element_filter(*type_it).getSize();
FEEngine & interface_engine =
this->MaterialReinforcement<dim>::model->getFEEngine("EmbeddedInterfaceFEEngine");
UInt nb_quad = interface_engine.getNbQuadraturePoints(*type_it);
this->ConstLaw::gradu.alloc(nb_elements * nb_quad, 1, *type_it, _not_ghost);
this->ConstLaw::stress.alloc(nb_elements * nb_quad, 1, *type_it, _not_ghost);
this->ConstLaw::delta_T.alloc(nb_elements * nb_quad, 1, *type_it, _not_ghost);
this->ConstLaw::sigma_th.alloc(nb_elements * nb_quad, 1, *type_it, _not_ghost);
this->ConstLaw::potential_energy.alloc(nb_elements * nb_quad, 1, *type_it, _not_ghost);
//this->ConstLaw::gradu.alloc(0, dim * dim, *type_it, _ghost);
}
}
/* -------------------------------------------------------------------------- */
template<UInt dim, class ConstLaw>
void MaterialReinforcementTemplate<dim, ConstLaw>::computeGradU(const ElementType & el_type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
MaterialReinforcement<dim>::computeGradU(el_type, ghost_type);
const UInt voigt_size = Material::getTangentStiffnessVoigtSize(dim);
Array<Real>::matrix_iterator full_gradu_it =
this->MaterialReinforcement<dim>::gradu(el_type, ghost_type).begin(dim, dim);
Array<Real>::matrix_iterator full_gradu_end =
this->MaterialReinforcement<dim>::gradu(el_type, ghost_type).end(dim, dim);
Array<Real>::scalar_iterator gradu_it =
this->ConstLaw::gradu(el_type, ghost_type).begin();
Array<Real>::matrix_iterator cosines_it =
this->directing_cosines(el_type, ghost_type).begin(voigt_size, voigt_size);
for (; full_gradu_it != full_gradu_end ; ++full_gradu_it,
++gradu_it,
++cosines_it) {
Matrix<Real> & full_gradu = *full_gradu_it;
Real & gradu = *gradu_it;
Matrix<Real> & C = *cosines_it;
computeInterfaceGradUOnQuad(full_gradu, gradu, C);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim, class ConstLaw>
void MaterialReinforcementTemplate<dim, ConstLaw>::computeInterfaceGradUOnQuad(const Matrix<Real> & full_gradu,
Real & gradu,
const Matrix<Real> & C) {
const UInt voigt_size = Material::getTangentStiffnessVoigtSize(dim);
Matrix<Real> epsilon(dim, dim);
Vector<Real> e_voigt(voigt_size);
Vector<Real> e_interface_voigt(voigt_size);
epsilon = 0.5 * (full_gradu + full_gradu.transpose());
MaterialReinforcement<dim>::strainTensorToVoigtVector(epsilon, e_voigt);
e_interface_voigt.mul<false>(C, e_voigt);
gradu = e_interface_voigt(0);
}
/* -------------------------------------------------------------------------- */
template<UInt dim, class ConstLaw>
void MaterialReinforcementTemplate<dim, ConstLaw>::computeTangentModuli(const ElementType & el_type,
Array<Real> & tangent,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(tangent.getNbComponent() == 1, "Reinforcements only work in 1D");
ConstLaw::computeTangentModuli(el_type, tangent, ghost_type);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim, class ConstLaw>
void MaterialReinforcementTemplate<dim, ConstLaw>::computeStress(ElementType type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
ConstLaw::computeStress(type, ghost_type);
Array<Real>::matrix_iterator full_sigma_it =
this->MaterialReinforcement<dim>::stress(type, ghost_type).begin(dim, dim);
Array<Real>::matrix_iterator full_sigma_end =
this->MaterialReinforcement<dim>::stress(type, ghost_type).end(dim, dim);
Array<Real>::scalar_iterator sigma_it =
this->ConstLaw::stress(type, ghost_type).begin();
+ Array<Real>::scalar_iterator pre_stress_it =
+ this->MaterialReinforcement<dim>::pre_stress(type, ghost_type).begin();
- for (; full_sigma_it != full_sigma_end ; ++full_sigma_it, ++sigma_it) {
+ for (; full_sigma_it != full_sigma_end ; ++full_sigma_it, ++sigma_it, ++pre_stress_it) {
Matrix<Real> & sigma = *full_sigma_it;
- sigma(0, 0) = *sigma_it;
+ sigma(0, 0) = *sigma_it + *pre_stress_it;
}
AKANTU_DEBUG_OUT();
}
diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc
index 78de0b7d5..d9dad7b21 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model.cc
@@ -1,1833 +1,1832 @@
/**
* @file solid_mechanics_model.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Tue Jul 27 2010
* @date last modification: Fri Sep 19 2014
*
* @brief Implementation of the SolidMechanicsModel class
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_math.hh"
#include "aka_common.hh"
#include "solid_mechanics_model.hh"
#include "group_manager_inline_impl.cc"
#include "dumpable_inline_impl.hh"
#include "integration_scheme_2nd_order.hh"
#include "element_group.hh"
#include "static_communicator.hh"
#include "dof_synchronizer.hh"
#include "element_group.hh"
#include <cmath>
#ifdef AKANTU_USE_MUMPS
#include "solver_mumps.hh"
#endif
#ifdef AKANTU_USE_PETSC
#include "solver_petsc.hh"
#include "petsc_matrix.hh"
#endif
#ifdef AKANTU_USE_IOHELPER
# include "dumper_field.hh"
# include "dumper_paraview.hh"
# include "dumper_homogenizing_field.hh"
# include "dumper_material_internal_field.hh"
# include "dumper_elemental_field.hh"
# include "dumper_material_padders.hh"
# include "dumper_element_partition.hh"
# include "dumper_iohelper.hh"
#endif
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
const SolidMechanicsModelOptions default_solid_mechanics_model_options(_explicit_lumped_mass, false);
/* -------------------------------------------------------------------------- */
/**
* 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
*/
SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh,
UInt dim,
const ID & id,
const MemoryID & memory_id) :
Model(mesh, dim, id, memory_id),
BoundaryCondition<SolidMechanicsModel>(),
time_step(NAN), f_m2a(1.0),
mass_matrix(NULL),
velocity_damping_matrix(NULL),
stiffness_matrix(NULL),
jacobian_matrix(NULL),
material_index("material index", id),
material_local_numbering("material local numbering", id),
material_selector(new DefaultMaterialSelector(material_index)),
is_default_material_selector(true),
integrator(NULL),
increment_flag(false), solver(NULL),
synch_parallel(NULL),
are_materials_instantiated(false) {
AKANTU_DEBUG_IN();
createSynchronizerRegistry(this);
registerFEEngineObject<MyFEEngineType>("SolidMechanicsFEEngine", mesh, spatial_dimension);
this->displacement = NULL;
this->mass = NULL;
this->velocity = NULL;
this->acceleration = NULL;
this->force = NULL;
this->residual = NULL;
this->blocked_dofs = NULL;
this->increment = NULL;
this->increment_acceleration = NULL;
this->dof_synchronizer = NULL;
this->previous_displacement = NULL;
materials.clear();
mesh.registerEventHandler(*this);
#if defined(AKANTU_USE_IOHELPER)
this->mesh.registerDumper<DumperParaview>("paraview_all", id, true);
this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular);
#endif
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
SolidMechanicsModel::~SolidMechanicsModel() {
AKANTU_DEBUG_IN();
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
delete *mat_it;
}
materials.clear();
delete integrator;
delete solver;
delete mass_matrix;
delete velocity_damping_matrix;
if(stiffness_matrix && stiffness_matrix != jacobian_matrix)
delete stiffness_matrix;
delete jacobian_matrix;
delete synch_parallel;
if(is_default_material_selector) {
delete material_selector;
material_selector = NULL;
}
AKANTU_DEBUG_OUT();
}
void SolidMechanicsModel::setTimeStep(Real time_step) {
this->time_step = time_step;
#if defined(AKANTU_USE_IOHELPER)
this->mesh.getDumper().setTimeStep(time_step);
#endif
}
/* -------------------------------------------------------------------------- */
/* Initialisation */
/* -------------------------------------------------------------------------- */
/**
* 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 material_file the file containing the materials to use
* @param method the analysis method wanted. See the akantu::AnalysisMethod for
* the different possibilities
*/
void SolidMechanicsModel::initFull(const ModelOptions & options) {
Model::initFull(options);
const SolidMechanicsModelOptions & smm_options =
dynamic_cast<const SolidMechanicsModelOptions &>(options);
method = smm_options.analysis_method;
// initialize the vectors
initArrays();
// set the initial condition to 0
force->clear();
velocity->clear();
acceleration->clear();
displacement->clear();
// initialize pcb
if(pbc_pair.size()!=0)
initPBC();
// initialize the time integration schemes
switch(method) {
case _explicit_lumped_mass:
initExplicit();
break;
case _explicit_consistent_mass:
initSolver();
initExplicit();
break;
case _implicit_dynamic:
initImplicit(true);
break;
case _static:
initImplicit(false);
break;
default:
AKANTU_EXCEPTION("analysis method not recognised by SolidMechanicsModel");
break;
}
// initialize the materials
if(this->parser->getLastParsedFile() != "") {
instantiateMaterials();
}
if(!smm_options.no_init_materials) {
initMaterials();
}
if(increment_flag)
initBC(*this, *displacement, *increment, *force);
else
initBC(*this, *displacement, *force);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initParallel(MeshPartition * partition,
DataAccessor * data_accessor) {
AKANTU_DEBUG_IN();
if (data_accessor == NULL) data_accessor = this;
synch_parallel = &createParallelSynch(partition,data_accessor);
synch_registry->registerSynchronizer(*synch_parallel, _gst_material_id);
synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_mass);
synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_stress);
synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_boundary);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initFEEngineBoundary() {
FEEngine & fem_boundary = getFEEngineBoundary();
fem_boundary.initShapeFunctions(_not_ghost);
fem_boundary.initShapeFunctions(_ghost);
fem_boundary.computeNormalsOnControlPoints(_not_ghost);
fem_boundary.computeNormalsOnControlPoints(_ghost);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initExplicit(AnalysisMethod analysis_method) {
AKANTU_DEBUG_IN();
//in case of switch from implicit to explicit
if(!this->isExplicit())
method = analysis_method;
if (integrator) delete integrator;
integrator = new CentralDifference();
UInt nb_nodes = acceleration->getSize();
UInt nb_degree_of_freedom = acceleration->getNbComponent();
std::stringstream sstr; sstr << id << ":increment_acceleration";
increment_acceleration = &(alloc<Real>(sstr.str(), nb_nodes, nb_degree_of_freedom, Real()));
AKANTU_DEBUG_OUT();
}
void SolidMechanicsModel::initArraysPreviousDisplacment() {
AKANTU_DEBUG_IN();
SolidMechanicsModel::setIncrementFlagOn();
UInt nb_nodes = mesh.getNbNodes();
std::stringstream sstr_disp_t;
sstr_disp_t << id << ":previous_displacement";
previous_displacement = &(alloc<Real > (sstr_disp_t.str(), nb_nodes, spatial_dimension, 0.));
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Allocate all the needed vectors. By default their are not necessarily set to
* 0
*
*/
void SolidMechanicsModel::initArrays() {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
std::stringstream sstr_disp; sstr_disp << id << ":displacement";
// std::stringstream sstr_mass; sstr_mass << id << ":mass";
std::stringstream sstr_velo; sstr_velo << id << ":velocity";
std::stringstream sstr_acce; sstr_acce << id << ":acceleration";
std::stringstream sstr_forc; sstr_forc << id << ":force";
std::stringstream sstr_resi; sstr_resi << id << ":residual";
std::stringstream sstr_boun; sstr_boun << id << ":blocked_dofs";
displacement = &(alloc<Real>(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
// mass = &(alloc<Real>(sstr_mass.str(), nb_nodes, spatial_dimension, 0));
velocity = &(alloc<Real>(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
acceleration = &(alloc<Real>(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
force = &(alloc<Real>(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
residual = &(alloc<Real>(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
blocked_dofs = &(alloc<bool>(sstr_boun.str(), nb_nodes, spatial_dimension, false));
std::stringstream sstr_curp; sstr_curp << id << ":current_position";
current_position = &(alloc<Real>(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE));
for(UInt g = _not_ghost; g <= _ghost; ++g) {
GhostType gt = (GhostType) g;
Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined);
Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined);
for(; it != end; ++it) {
UInt nb_element = mesh.getNbElement(*it, gt);
material_index.alloc(nb_element, 1, *it, gt);
material_local_numbering.alloc(nb_element, 1, *it, gt);
}
}
dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension);
dof_synchronizer->initLocalDOFEquationNumbers();
dof_synchronizer->initGlobalDOFEquationNumbers();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* 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::initPBC() {
Model::initPBC();
registerPBCSynchronizer();
// as long as there are ones on the diagonal of the matrix, we can put boudandary true for slaves
std::map<UInt, UInt>::iterator it = pbc_pair.begin();
std::map<UInt, UInt>::iterator end = pbc_pair.end();
UInt dim = mesh.getSpatialDimension();
while(it != end) {
for (UInt i=0; i<dim; ++i)
(*blocked_dofs)((*it).first,i) = true;
++it;
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::registerPBCSynchronizer(){
PBCSynchronizer * synch = new PBCSynchronizer(pbc_pair);
synch_registry->registerSynchronizer(*synch, _gst_smm_uv);
synch_registry->registerSynchronizer(*synch, _gst_smm_mass);
synch_registry->registerSynchronizer(*synch, _gst_smm_res);
synch_registry->registerSynchronizer(*synch, _gst_for_dump);
changeLocalEquationNumberForPBC(pbc_pair, mesh.getSpatialDimension());
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateCurrentPosition() {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
current_position->resize(nb_nodes);
Real * current_position_val = current_position->storage();
Real * position_val = mesh.getNodes().storage();
Real * displacement_val = displacement->storage();
/// compute current_position = initial_position + displacement
memcpy(current_position_val, position_val, nb_nodes*spatial_dimension*sizeof(Real));
for (UInt n = 0; n < nb_nodes*spatial_dimension; ++n) {
*current_position_val++ += *displacement_val++;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initializeUpdateResidualData() {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
residual->resize(nb_nodes);
/// copy the forces in residual for boundary conditions
memcpy(residual->storage(), force->storage(), nb_nodes*spatial_dimension*sizeof(Real));
// start synchronization
synch_registry->asynchronousSynchronize(_gst_smm_uv);
synch_registry->waitEndSynchronize(_gst_smm_uv);
updateCurrentPosition();
AKANTU_DEBUG_OUT();
}
/*----------------------------------------------------------------------------*/
void SolidMechanicsModel::reInitialize()
{
}
/* -------------------------------------------------------------------------- */
/* Explicit scheme */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/**
* This function compute the second member of the motion equation. That is to
* say the sum of forces @f$ r = F_{ext} - F_{int} @f$. @f$ F_{ext} @f$ is given
* by the user in the force vector, and @f$ F_{int} @f$ is computed as @f$
* F_{int} = \int_{\Omega} N \sigma d\Omega@f$
*
*/
void SolidMechanicsModel::updateResidual(bool need_initialize) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the internal forces");
// f = f_ext - f_int
// f = f_ext
if(need_initialize) initializeUpdateResidualData();
AKANTU_DEBUG_INFO("Compute local stresses");
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllStresses(_not_ghost);
}
#ifdef AKANTU_DAMAGE_NON_LOCAL
/* ------------------------------------------------------------------------ */
/* Computation of the non local part */
synch_registry->asynchronousSynchronize(_gst_mnl_for_average);
AKANTU_DEBUG_INFO("Compute non local stresses for local elements");
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllNonLocalStresses(_not_ghost);
}
AKANTU_DEBUG_INFO("Wait distant non local stresses");
synch_registry->waitEndSynchronize(_gst_mnl_for_average);
AKANTU_DEBUG_INFO("Compute non local stresses for ghosts elements");
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllNonLocalStresses(_ghost);
}
#endif
/* ------------------------------------------------------------------------ */
/* assembling the forces internal */
// communicate the stress
AKANTU_DEBUG_INFO("Send data for residual assembly");
synch_registry->asynchronousSynchronize(_gst_smm_stress);
AKANTU_DEBUG_INFO("Assemble residual for local elements");
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.assembleResidual(_not_ghost);
}
AKANTU_DEBUG_INFO("Wait distant stresses");
// finalize communications
synch_registry->waitEndSynchronize(_gst_smm_stress);
AKANTU_DEBUG_INFO("Assemble residual for ghost elements");
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.assembleResidual(_ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::computeStresses() {
if (isExplicit()) {
// start synchronization
synch_registry->asynchronousSynchronize(_gst_smm_uv);
synch_registry->waitEndSynchronize(_gst_smm_uv);
// compute stresses on all local elements for each materials
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllStresses(_not_ghost);
}
/* ------------------------------------------------------------------------ */
#ifdef AKANTU_DAMAGE_NON_LOCAL
/* Computation of the non local part */
synch_registry->asynchronousSynchronize(_gst_mnl_for_average);
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllNonLocalStresses(_not_ghost);
}
synch_registry->waitEndSynchronize(_gst_mnl_for_average);
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllNonLocalStresses(_ghost);
}
#endif
} else {
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllStressesFromTangentModuli(_not_ghost);
}
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateResidualInternal() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Update the residual");
// f = f_ext - f_int - Ma - Cv = r - Ma - Cv;
if(method != _static) {
// f -= Ma
if(mass_matrix) {
// if full mass_matrix
Array<Real> * Ma = new Array<Real>(*acceleration, true, "Ma");
*Ma *= *mass_matrix;
/// \todo check unit conversion for implicit dynamics
// *Ma /= f_m2a
*residual -= *Ma;
delete Ma;
} else if (mass) {
// else lumped mass
UInt nb_nodes = acceleration->getSize();
UInt nb_degree_of_freedom = acceleration->getNbComponent();
Real * mass_val = mass->storage();
Real * accel_val = acceleration->storage();
Real * res_val = residual->storage();
bool * blocked_dofs_val = blocked_dofs->storage();
for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) {
if(!(*blocked_dofs_val)) {
*res_val -= *accel_val * *mass_val /f_m2a;
}
blocked_dofs_val++;
res_val++;
mass_val++;
accel_val++;
}
} else {
AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix.");
}
// f -= Cv
if(velocity_damping_matrix) {
Array<Real> * Cv = new Array<Real>(*velocity);
*Cv *= *velocity_damping_matrix;
*residual -= *Cv;
delete Cv;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateAcceleration() {
AKANTU_DEBUG_IN();
updateResidualInternal();
if(method == _explicit_lumped_mass) {
/* residual = residual_{n+1} - M * acceleration_n therefore
solution = increment acceleration not acceleration */
solveLumped(*increment_acceleration,
*mass,
*residual,
*blocked_dofs,
f_m2a);
} else if (method == _explicit_consistent_mass) {
solve<NewmarkBeta::_acceleration_corrector>(*increment_acceleration);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::solveLumped(Array<Real> & x,
const Array<Real> & A,
const Array<Real> & b,
const Array<bool> & blocked_dofs,
Real alpha) {
Real * A_val = A.storage();
Real * b_val = b.storage();
Real * x_val = x.storage();
bool * blocked_dofs_val = blocked_dofs.storage();
UInt nb_degrees_of_freedom = x.getSize() * x.getNbComponent();
for (UInt n = 0; n < nb_degrees_of_freedom; ++n) {
if(!(*blocked_dofs_val)) {
*x_val = alpha * (*b_val / *A_val);
}
x_val++;
A_val++;
b_val++;
blocked_dofs_val++;
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::explicitPred() {
AKANTU_DEBUG_IN();
if(increment_flag) {
if(previous_displacement) increment->copy(*previous_displacement);
else increment->copy(*displacement);
}
AKANTU_DEBUG_ASSERT(integrator,"itegrator should have been allocated: "
<< "have called initExplicit ? "
<< "or initImplicit ?");
integrator->integrationSchemePred(time_step,
*displacement,
*velocity,
*acceleration,
*blocked_dofs);
if(increment_flag) {
Real * inc_val = increment->storage();
Real * dis_val = displacement->storage();
UInt nb_degree_of_freedom = displacement->getSize() * displacement->getNbComponent();
for (UInt n = 0; n < nb_degree_of_freedom; ++n) {
*inc_val = *dis_val - *inc_val;
inc_val++;
dis_val++;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::explicitCorr() {
AKANTU_DEBUG_IN();
integrator->integrationSchemeCorrAccel(time_step,
*displacement,
*velocity,
*acceleration,
*blocked_dofs,
*increment_acceleration);
if(previous_displacement) previous_displacement->copy(*displacement);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::solveStep() {
AKANTU_DEBUG_IN();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeSolveStepEvent(method));
this->explicitPred();
this->updateResidual();
this->updateAcceleration();
this->explicitCorr();
EventManager::sendEvent(SolidMechanicsModelEvent::AfterSolveStepEvent(method));
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/* Implicit scheme */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/**
* Initialize the solver and create the sparse matrices needed.
*
*/
void SolidMechanicsModel::initSolver(__attribute__((unused)) SolverOptions & options) {
#if !defined(AKANTU_USE_MUMPS) && !defined(AKANTU_USE_PETSC)// or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake
AKANTU_DEBUG_ERROR("You should at least activate one solver.");
#else
UInt nb_global_nodes = mesh.getNbGlobalNodes();
delete jacobian_matrix;
std::stringstream sstr; sstr << id << ":jacobian_matrix";
#ifdef AKANTU_USE_PETSC
jacobian_matrix = new PETScMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id);
#else
jacobian_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id);
#endif //AKANTU_USE PETSC
jacobian_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension);
if (!isExplicit()) {
delete stiffness_matrix;
std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix";
#ifdef AKANTU_USE_PETSC
stiffness_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id);
stiffness_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension);
#else
stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id);
#endif //AKANTU_USE_PETSC
}
delete solver;
std::stringstream sstr_solv; sstr_solv << id << ":solver";
#ifdef AKANTU_USE_PETSC
solver = new SolverPETSc(*jacobian_matrix, sstr_solv.str());
#elif defined(AKANTU_USE_MUMPS)
solver = new SolverMumps(*jacobian_matrix, sstr_solv.str());
dof_synchronizer->initScatterGatherCommunicationScheme();
#else
AKANTU_DEBUG_ERROR("You should at least activate one solver.");
#endif //AKANTU_USE_MUMPS
if(solver)
solver->initialize(options);
#endif //AKANTU_HAS_SOLVER
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initJacobianMatrix() {
#if defined(AKANTU_USE_MUMPS) && !defined(AKANTU_USE_PETSC)
// @todo make it more flexible: this is an ugly patch to treat the case of non
// fix profile of the K matrix
delete jacobian_matrix;
std::stringstream sstr_sti; sstr_sti << id << ":jacobian_matrix";
jacobian_matrix = new SparseMatrix(*stiffness_matrix, sstr_sti.str(), memory_id);
std::stringstream sstr_solv; sstr_solv << id << ":solver";
delete solver;
solver = new SolverMumps(*jacobian_matrix, sstr_solv.str());
if(solver)
solver->initialize(_solver_no_options);
#else
AKANTU_DEBUG_ERROR("You need to activate the solver mumps.");
#endif
}
/* -------------------------------------------------------------------------- */
/**
* Initialize the implicit solver, either for dynamic or static cases,
*
* @param dynamic
*/
void SolidMechanicsModel::initImplicit(bool dynamic, SolverOptions & solver_options) {
AKANTU_DEBUG_IN();
method = dynamic ? _implicit_dynamic : _static;
if (!increment) setIncrementFlagOn();
initSolver(solver_options);
if(method == _implicit_dynamic) {
if(integrator) delete integrator;
integrator = new TrapezoidalRule2();
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initialAcceleration() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Solving Ma = f");
Solver * acc_solver = NULL;
std::stringstream sstr; sstr << id << ":tmp_mass_matrix";
SparseMatrix * tmp_mass = new SparseMatrix(*mass_matrix, sstr.str(), memory_id);
#ifdef AKANTU_USE_MUMPS
std::stringstream sstr_solver; sstr << id << ":solver_mass_matrix";
acc_solver = new SolverMumps(*mass_matrix, sstr_solver.str());
dof_synchronizer->initScatterGatherCommunicationScheme();
#else
AKANTU_DEBUG_ERROR("You should at least activate one solver.");
#endif //AKANTU_USE_MUMPS
acc_solver->initialize();
tmp_mass->applyBoundary(*blocked_dofs);
acc_solver->setRHS(*residual);
acc_solver->solve(*acceleration);
delete acc_solver;
delete tmp_mass;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleStiffnessMatrix() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the new stiffness matrix.");
stiffness_matrix->clear();
// call compute stiffness matrix on each local elements
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->assembleStiffnessMatrix(_not_ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
SparseMatrix & SolidMechanicsModel::initVelocityDampingMatrix() {
if(!velocity_damping_matrix)
velocity_damping_matrix =
new SparseMatrix(*jacobian_matrix, id + ":velocity_damping_matrix", memory_id);
return *velocity_damping_matrix;
}
/* -------------------------------------------------------------------------- */
template<>
bool SolidMechanicsModel::testConvergence<_scc_increment>(Real tolerance, Real & error){
AKANTU_DEBUG_IN();
UInt nb_nodes = displacement->getSize();
UInt nb_degree_of_freedom = displacement->getNbComponent();
error = 0;
Real norm[2] = {0., 0.};
Real * increment_val = increment->storage();
bool * blocked_dofs_val = blocked_dofs->storage();
Real * displacement_val = displacement->storage();
for (UInt n = 0; n < nb_nodes; ++n) {
bool is_local_node = mesh.isLocalOrMasterNode(n);
for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
if(!(*blocked_dofs_val) && is_local_node) {
norm[0] += *increment_val * *increment_val;
norm[1] += *displacement_val * *displacement_val;
}
blocked_dofs_val++;
increment_val++;
displacement_val++;
}
}
StaticCommunicator::getStaticCommunicator().allReduce(norm, 2, _so_sum);
norm[0] = sqrt(norm[0]);
norm[1] = sqrt(norm[1]);
AKANTU_DEBUG_ASSERT(!Math::isnan(norm[0]), "Something goes wrong in the solve phase");
if (norm[1] < Math::getTolerance()) {
error = norm[0];
AKANTU_DEBUG_OUT();
// cout<<"Error 1: "<<error<<endl;
return error < tolerance;
}
AKANTU_DEBUG_OUT();
if(norm[1] > Math::getTolerance())
error = norm[0] / norm[1];
else
error = norm[0]; //In case the total displacement is zero!
// cout<<"Error 2: "<<error<<endl;
return (error < tolerance);
}
/* -------------------------------------------------------------------------- */
template<>
bool SolidMechanicsModel::testConvergence<_scc_residual>(Real tolerance, Real & norm) {
AKANTU_DEBUG_IN();
-
-
UInt nb_nodes = residual->getSize();
+ UInt nb_degree_of_freedom = displacement->getNbComponent();
norm = 0;
Real * residual_val = residual->storage();
bool * blocked_dofs_val = blocked_dofs->storage();
for (UInt n = 0; n < nb_nodes; ++n) {
bool is_local_node = mesh.isLocalOrMasterNode(n);
if(is_local_node) {
- for (UInt d = 0; d < spatial_dimension; ++d) {
+ for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
if(!(*blocked_dofs_val)) {
norm += *residual_val * *residual_val;
}
blocked_dofs_val++;
residual_val++;
}
} else {
blocked_dofs_val += spatial_dimension;
residual_val += spatial_dimension;
}
}
StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum);
norm = sqrt(norm);
AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase");
AKANTU_DEBUG_OUT();
return (norm < tolerance);
}
/* -------------------------------------------------------------------------- */
template<>
bool SolidMechanicsModel::testConvergence<_scc_residual_mass_wgh>(Real tolerance,
Real & norm) {
AKANTU_DEBUG_IN();
UInt nb_nodes = residual->getSize();
norm = 0;
Real * residual_val = residual->storage();
Real * mass_val = this->mass->storage();
bool * blocked_dofs_val = blocked_dofs->storage();
for (UInt n = 0; n < nb_nodes; ++n) {
bool is_local_node = mesh.isLocalOrMasterNode(n);
if(is_local_node) {
for (UInt d = 0; d < spatial_dimension; ++d) {
if(!(*blocked_dofs_val)) {
norm += *residual_val * *residual_val/(*mass_val * *mass_val);
}
blocked_dofs_val++;
residual_val++;
mass_val++;
}
} else {
blocked_dofs_val += spatial_dimension;
residual_val += spatial_dimension;
mass_val += spatial_dimension;
}
}
StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum);
norm = sqrt(norm);
AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase");
AKANTU_DEBUG_OUT();
return (norm < tolerance);
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::testConvergenceResidual(Real tolerance){
AKANTU_DEBUG_IN();
Real error=0;
bool res = this->testConvergence<_scc_residual>(tolerance, error);
AKANTU_DEBUG_OUT();
return res;
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::testConvergenceResidual(Real tolerance, Real & error){
AKANTU_DEBUG_IN();
bool res = this->testConvergence<_scc_residual>(tolerance, error);
AKANTU_DEBUG_OUT();
return res;
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance){
AKANTU_DEBUG_IN();
Real error=0;
bool res = this->testConvergence<_scc_increment>(tolerance, error);
AKANTU_DEBUG_OUT();
return res;
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance, Real & error){
AKANTU_DEBUG_IN();
bool res = this->testConvergence<_scc_increment>(tolerance, error);
AKANTU_DEBUG_OUT();
return res;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::implicitPred() {
AKANTU_DEBUG_IN();
if(previous_displacement) previous_displacement->copy(*displacement);
if(method == _implicit_dynamic)
integrator->integrationSchemePred(time_step,
*displacement,
*velocity,
*acceleration,
*blocked_dofs);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::implicitCorr() {
AKANTU_DEBUG_IN();
if(method == _implicit_dynamic) {
integrator->integrationSchemeCorrDispl(time_step,
*displacement,
*velocity,
*acceleration,
*blocked_dofs,
*increment);
} else {
UInt nb_nodes = displacement->getSize();
UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes;
Real * incr_val = increment->storage();
Real * disp_val = displacement->storage();
bool * boun_val = blocked_dofs->storage();
for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, ++boun_val){
*incr_val *= (1. - *boun_val);
*disp_val += *incr_val;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateIncrement() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(previous_displacement,"The previous displacement has to be initialized."
<< " Are you working with Finite or Ineslactic deformations?");
UInt nb_nodes = displacement->getSize();
UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes;
Real * incr_val = increment->storage();
Real * disp_val = displacement->storage();
Real * prev_disp_val = previous_displacement->storage();
for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, ++prev_disp_val)
*incr_val = (*disp_val - *prev_disp_val);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updatePreviousDisplacement() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(previous_displacement,"The previous displacement has to be initialized."
<< " Are you working with Finite or Ineslactic deformations?");
previous_displacement->copy(*displacement);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/* Information */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::synchronizeBoundaries() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized."
<< " Did you call initParallel?");
synch_registry->synchronize(_gst_smm_boundary);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::synchronizeResidual() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized."
<< " Did you call initPBC?");
synch_registry->synchronize(_gst_smm_res);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::setIncrementFlagOn() {
AKANTU_DEBUG_IN();
if(!increment) {
UInt nb_nodes = mesh.getNbNodes();
std::stringstream sstr_inc; sstr_inc << id << ":increment";
increment = &(alloc<Real>(sstr_inc.str(), nb_nodes, spatial_dimension, 0.));
}
increment_flag = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getStableTimeStep() {
AKANTU_DEBUG_IN();
Real min_dt = getStableTimeStep(_not_ghost);
/// reduction min over all processors
StaticCommunicator::getStaticCommunicator().allReduce(&min_dt, 1, _so_min);
AKANTU_DEBUG_OUT();
return min_dt;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) {
AKANTU_DEBUG_IN();
Material ** mat_val = &(materials.at(0));
Real min_dt = std::numeric_limits<Real>::max();
updateCurrentPosition();
Element elem;
elem.ghost_type = ghost_type;
elem.kind = _ek_regular;
Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type);
for(; it != end; ++it) {
elem.type = *it;
UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it);
UInt nb_element = mesh.getNbElement(*it);
Array<UInt>::const_scalar_iterator mat_indexes = material_index(*it, ghost_type).begin();
Array<UInt>::const_scalar_iterator mat_loc_num = material_local_numbering(*it, ghost_type).begin();
Array<Real> X(0, nb_nodes_per_element*spatial_dimension);
FEEngine::extractNodalToElementField(mesh, *current_position,
X, *it, _not_ghost);
Array<Real>::matrix_iterator X_el =
X.begin(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, *it);
Real el_c = mat_val[*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::getPotentialEnergy() {
AKANTU_DEBUG_IN();
Real energy = 0.;
/// call update residual on each local elements
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
energy += (*mat_it)->getPotentialEnergy();
}
/// reduction sum over all processors
StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum);
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getKineticEnergy() {
AKANTU_DEBUG_IN();
if (!mass && !mass_matrix)
AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix.");
Real ekin = 0.;
UInt nb_nodes = mesh.getNbNodes();
Real * vel_val = velocity->storage();
Real * mass_val = mass->storage();
for (UInt n = 0; n < nb_nodes; ++n) {
Real mv2 = 0;
bool is_local_node = mesh.isLocalOrMasterNode(n);
bool is_not_pbc_slave_node = !isPBCSlaveNode(n);
bool count_node = is_local_node && is_not_pbc_slave_node;
for (UInt i = 0; i < spatial_dimension; ++i) {
if (count_node)
mv2 += *vel_val * *vel_val * *mass_val;
vel_val++;
mass_val++;
}
ekin += mv2;
}
StaticCommunicator::getStaticCommunicator().allReduce(&ekin, 1, _so_sum);
AKANTU_DEBUG_OUT();
return ekin * .5;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getKineticEnergy(const ElementType & type, UInt index) {
AKANTU_DEBUG_IN();
UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(type);
Array<Real> vel_on_quad(nb_quadrature_points, spatial_dimension);
Array<UInt> filter_element(1, 1, index);
getFEEngine().interpolateOnQuadraturePoints(*velocity, vel_on_quad,
spatial_dimension,
type, _not_ghost,
filter_element);
Array<Real>::vector_iterator vit = vel_on_quad.begin(spatial_dimension);
Array<Real>::vector_iterator vend = vel_on_quad.end(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();
Real * velo = velocity->storage();
Real * forc = force->storage();
Real * resi = residual->storage();
bool * boun = blocked_dofs->storage();
Real work = 0.;
UInt nb_nodes = mesh.getNbNodes();
for (UInt n = 0; n < nb_nodes; ++n) {
bool is_local_node = mesh.isLocalOrMasterNode(n);
bool is_not_pbc_slave_node = !isPBCSlaveNode(n);
bool count_node = is_local_node && is_not_pbc_slave_node;
for (UInt i = 0; i < spatial_dimension; ++i) {
if (count_node) {
if(*boun)
work -= *resi * *velo * time_step;
else
work += *forc * *velo * time_step;
}
++velo;
++forc;
++resi;
++boun;
}
}
StaticCommunicator::getStaticCommunicator().allReduce(&work, 1, _so_sum);
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.;
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
energy += (*mat_it)->getEnergy(energy_id);
}
/// reduction sum over all processors
StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_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);
}
std::vector<Material *>::iterator mat_it;
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->getFEEngine().initShapeFunctions(_not_ghost);
this->getFEEngine().initShapeFunctions(_ghost);
for(UInt g = _not_ghost; g <= _ghost; ++g) {
GhostType gt = (GhostType) g;
Mesh::type_iterator it = this->mesh.firstType(spatial_dimension, gt, _ek_not_defined);
Mesh::type_iterator end = this->mesh.lastType(spatial_dimension, gt, _ek_not_defined);
for(; it != end; ++it) {
UInt nb_element = this->mesh.getNbElement(*it, gt);
if(!material_index.exists(*it, gt)) {
this->material_index .alloc(nb_element, 1, *it, gt);
this->material_local_numbering.alloc(nb_element, 1, *it, gt);
} else {
this->material_index (*it, gt).resize(nb_element);
this->material_local_numbering(*it, gt).resize(nb_element);
}
}
}
Array<Element>::const_iterator<Element> it = element_list.begin();
Array<Element>::const_iterator<Element> end = element_list.end();
ElementTypeMapArray<UInt> filter("new_element_filter", this->getID());
for (UInt el = 0; it != end; ++it, ++el) {
const Element & elem = *it;
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->assignMaterialToElements(&filter);
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->onElementsAdded(element_list, event);
}
if(method == _explicit_lumped_mass) this->assembleMassLumped();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onElementsRemoved(__attribute__((unused)) const Array<Element> & element_list,
const ElementTypeMapArray<UInt> & new_numbering,
const RemovedElementsEvent & event) {
this->getFEEngine().initShapeFunctions(_not_ghost);
this->getFEEngine().initShapeFunctions(_ghost);
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->onElementsRemoved(element_list, new_numbering, event);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onNodesAdded(const Array<UInt> & nodes_list,
__attribute__((unused)) const NewNodesEvent & event) {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
if(displacement) displacement->resize(nb_nodes);
if(mass ) mass ->resize(nb_nodes);
if(velocity ) velocity ->resize(nb_nodes);
if(acceleration) acceleration->resize(nb_nodes);
if(force ) force ->resize(nb_nodes);
if(residual ) residual ->resize(nb_nodes);
if(blocked_dofs) blocked_dofs->resize(nb_nodes);
if(previous_displacement) previous_displacement->resize(nb_nodes);
if(increment_acceleration) increment_acceleration->resize(nb_nodes);
if(increment) increment->resize(nb_nodes);
if(current_position) current_position->resize(nb_nodes);
delete dof_synchronizer;
dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension);
dof_synchronizer->initLocalDOFEquationNumbers();
dof_synchronizer->initGlobalDOFEquationNumbers();
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->onNodesAdded(nodes_list, event);
}
if (method != _explicit_lumped_mass) {
this->initSolver();
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onNodesRemoved(__attribute__((unused)) const Array<UInt> & element_list,
const Array<UInt> & new_numbering,
__attribute__((unused)) const RemovedNodesEvent & event) {
if(displacement) mesh.removeNodesFromArray(*displacement, new_numbering);
if(mass ) mesh.removeNodesFromArray(*mass , new_numbering);
if(velocity ) mesh.removeNodesFromArray(*velocity , new_numbering);
if(acceleration) mesh.removeNodesFromArray(*acceleration, new_numbering);
if(force ) mesh.removeNodesFromArray(*force , new_numbering);
if(residual ) mesh.removeNodesFromArray(*residual , new_numbering);
if(blocked_dofs) mesh.removeNodesFromArray(*blocked_dofs, new_numbering);
if(increment_acceleration) mesh.removeNodesFromArray(*increment_acceleration, new_numbering);
if(increment) mesh.removeNodesFromArray(*increment , new_numbering);
delete dof_synchronizer;
dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension);
dof_synchronizer->initLocalDOFEquationNumbers();
dof_synchronizer->initGlobalDOFEquationNumbers();
if (method != _explicit_lumped_mass) {
this->initSolver();
}
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::isInternal(const std::string & field_name,
const ElementKind & element_kind){
bool is_internal = false;
/// check if at least one material contains field_id as an internal
for (UInt m = 0; m < materials.size() && !is_internal; ++m) {
is_internal |= materials[m]->isInternal(field_name, element_kind);
}
return is_internal;
}
/* -------------------------------------------------------------------------- */
ElementTypeMap<UInt> SolidMechanicsModel::getInternalDataPerElem(const std::string & field_name,
const ElementKind & element_kind){
if (!(this->isInternal(field_name,element_kind))) AKANTU_EXCEPTION("unknown internal " << field_name);
for (UInt m = 0; m < materials.size() ; ++m) {
if (materials[m]->isInternal(field_name, element_kind))
return materials[m]->getInternalDataPerElem(field_name,element_kind);
}
return ElementTypeMap<UInt>();
}
/* -------------------------------------------------------------------------- */
ElementTypeMapArray<Real> & SolidMechanicsModel::flattenInternal(const std::string & field_name,
const ElementKind & kind, const GhostType ghost_type){
std::pair<std::string,ElementKind> key(field_name,kind);
if (this->registered_internals.count(key) == 0){
this->registered_internals[key] =
new ElementTypeMapArray<Real>(field_name,this->id);
}
ElementTypeMapArray<Real> * internal_flat = this->registered_internals[key];
for (UInt m = 0; m < materials.size(); ++m)
materials[m]->flattenInternal(field_name,*internal_flat,ghost_type,kind);
return *internal_flat;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::flattenAllRegisteredInternals(const ElementKind & kind){
std::map<std::pair<std::string,ElementKind>,ElementTypeMapArray<Real> *> ::iterator it = this->registered_internals.begin();
std::map<std::pair<std::string,ElementKind>,ElementTypeMapArray<Real> *>::iterator end = this->registered_internals.end();
while (it != end){
if (kind == it->first.second)
this->flattenInternal(it->first.first,kind);
++it;
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onDump(){
this->flattenAllRegisteredInternals(_ek_regular);
}
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
dumper::Field * SolidMechanicsModel
::createElementalField(const std::string & field_name,
const std::string & group_name,
bool padding_flag,
const ElementKind & kind){
dumper::Field * field = NULL;
if(field_name == "partitions")
field = mesh.createElementalField<UInt, dumper::ElementPartitionField>(mesh.getConnectivities(),group_name,this->spatial_dimension,kind);
else if(field_name == "material_index")
field = mesh.createElementalField<UInt, Vector, dumper::ElementalField >(material_index,group_name,this->spatial_dimension,kind);
else {
// this copy of field_name is used to compute derivated data such as
// strain and von mises stress that are based on grad_u and stress
std::string field_name_copy(field_name);
if (field_name == "strain"
|| field_name == "Green strain"
|| field_name == "principal strain"
|| field_name == "principal Green strain")
field_name_copy = "grad_u";
else if (field_name == "Von Mises stress")
field_name_copy = "stress";
bool is_internal = this->isInternal(field_name_copy,kind);
if (is_internal) {
ElementTypeMap<UInt> nb_data_per_elem = this->getInternalDataPerElem(field_name_copy,kind);
ElementTypeMapArray<Real> & internal_flat = this->flattenInternal(field_name_copy,kind);
field = mesh.createElementalField<Real, dumper::InternalMaterialField>(internal_flat,
group_name,
this->spatial_dimension,kind,nb_data_per_elem);
if (field_name == "strain"){
dumper::ComputeStrain<false> * foo = new dumper::ComputeStrain<false>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field,*foo);
} else if (field_name == "Von Mises stress") {
dumper::ComputeVonMisesStress * foo = new dumper::ComputeVonMisesStress(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field,*foo);
} else if (field_name == "Green strain") {
dumper::ComputeStrain<true> * foo = new dumper::ComputeStrain<true>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field,*foo);
} else if (field_name == "principal strain") {
dumper::ComputePrincipalStrain<false> * foo = new dumper::ComputePrincipalStrain<false>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field,*foo);
} else if (field_name == "principal Green strain") {
dumper::ComputePrincipalStrain<true> * foo = new dumper::ComputePrincipalStrain<true>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field,*foo);
}
//treat the paddings
if (padding_flag){
if (field_name == "stress"){
if (this->spatial_dimension == 2) {
dumper::StressPadder<2> * foo = new dumper::StressPadder<2>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field,*foo);
}
} else if (field_name == "strain" || field_name == "Green strain"){
if (this->spatial_dimension == 2) {
dumper::StrainPadder<2> * foo = new dumper::StrainPadder<2>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field,*foo);
}
}
}
// homogenize the field
dumper::ComputeFunctorInterface * foo =
dumper::HomogenizerProxy::createHomogenizer(*field);
field = dumper::FieldComputeProxy::createFieldCompute(field,*foo);
}
}
return field;
}
/* -------------------------------------------------------------------------- */
dumper::Field * SolidMechanicsModel::createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
std::map<std::string,Array<Real>* > real_nodal_fields;
real_nodal_fields["displacement" ] = displacement;
real_nodal_fields["mass" ] = mass;
real_nodal_fields["velocity" ] = velocity;
real_nodal_fields["acceleration" ] = acceleration;
real_nodal_fields["force" ] = force;
real_nodal_fields["residual" ] = residual;
real_nodal_fields["increment" ] = increment;
dumper::Field * field = NULL;
if (padding_flag)
field = mesh.createNodalField(real_nodal_fields[field_name],group_name, 3);
else
field = mesh.createNodalField(real_nodal_fields[field_name],group_name);
return field;
}
/* -------------------------------------------------------------------------- */
dumper::Field * SolidMechanicsModel::createNodalFieldBool(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
std::map<std::string,Array<bool>* > uint_nodal_fields;
uint_nodal_fields["blocked_dofs" ] = blocked_dofs;
dumper::Field * field = NULL;
field = mesh.createNodalField(uint_nodal_fields[field_name],group_name);
return field;
}
/* -------------------------------------------------------------------------- */
#else
/* -------------------------------------------------------------------------- */
dumper::Field * SolidMechanicsModel
::createElementalField(const std::string & field_name,
const std::string & group_name,
bool padding_flag,
const ElementKind & kind){
return NULL;
}
/* -------------------------------------------------------------------------- */
dumper::Field * SolidMechanicsModel::createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
return NULL;
}
/* -------------------------------------------------------------------------- */
dumper::Field * SolidMechanicsModel::createNodalFieldBool(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
return NULL;
}
#endif
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(const std::string & dumper_name) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
synch_registry->synchronize(_gst_for_dump);
mesh.dump(dumper_name);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(const std::string & dumper_name, UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
synch_registry->synchronize(_gst_for_dump);
mesh.dump(dumper_name, step);
}
/* ------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(const std::string & dumper_name, Real time, UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
synch_registry->synchronize(_gst_for_dump);
mesh.dump(dumper_name, time, step);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump() {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(step);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(Real time, UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(time, step);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::computeCauchyStresses() {
AKANTU_DEBUG_IN();
// call compute stiffness matrix on each local elements
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
if(mat.isFiniteDeformation())
mat.computeAllCauchyStresses(_not_ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::saveStressAndStrainBeforeDamage() {
EventManager::sendEvent(SolidMechanicsModelEvent::BeginningOfDamageIterationEvent());
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateEnergiesAfterDamage() {
EventManager::sendEvent(SolidMechanicsModelEvent::AfterDamageEvent());
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::printself(std::ostream & stream, int indent) const {
std::string space;
for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
stream << space << "Solid Mechanics Model [" << std::endl;
stream << space << " + id : " << id << std::endl;
stream << space << " + spatial dimension : " << spatial_dimension << std::endl;
stream << space << " + fem [" << std::endl;
getFEEngine().printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + nodals information [" << std::endl;
displacement->printself(stream, indent + 2);
mass ->printself(stream, indent + 2);
velocity ->printself(stream, indent + 2);
acceleration->printself(stream, indent + 2);
force ->printself(stream, indent + 2);
residual ->printself(stream, indent + 2);
blocked_dofs->printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + material information [" << std::endl;
material_index.printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + materials [" << std::endl;
std::vector<Material *>::const_iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
const Material & mat = *(*mat_it);
mat.printself(stream, indent + 1);
}
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__
diff --git a/src/model/structural_mechanics/structural_mechanics_model.cc b/src/model/structural_mechanics/structural_mechanics_model.cc
index e4a0111be..2c41ab2ac 100644
--- a/src/model/structural_mechanics/structural_mechanics_model.cc
+++ b/src/model/structural_mechanics/structural_mechanics_model.cc
@@ -1,1147 +1,1147 @@
/**
* @file structural_mechanics_model.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Damien Spielmann <damien.spielmann@epfl.ch>
* @author Sébastien Hartmann <sebastien.hartmann@epfl.ch>
* @author Fabian Barras <fabian.barras@epfl.ch>
*
* @date creation: Fri Jul 15 2011
* @date last modification: Mon Sep 15 2014
*
* @brief Model implementation for StucturalMechanics elements
*
* @section LICENSE
*
* Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "structural_mechanics_model.hh"
#include "group_manager_inline_impl.cc"
#include "dumpable_inline_impl.hh"
#include "aka_math.hh"
#include "integration_scheme_2nd_order.hh"
#include "static_communicator.hh"
#include "sparse_matrix.hh"
#include "solver.hh"
#ifdef AKANTU_USE_MUMPS
#include "solver_mumps.hh"
#endif
#ifdef AKANTU_USE_IOHELPER
# include "dumper_paraview.hh"
# include "dumper_elemental_field.hh"
#endif
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
const StructuralMechanicsModelOptions default_structural_mechanics_model_options(_static);
/* -------------------------------------------------------------------------- */
StructuralMechanicsModel::StructuralMechanicsModel(Mesh & mesh,
UInt dim,
const ID & id,
const MemoryID & memory_id) :
Model(mesh, dim, id, memory_id),
time_step(NAN), f_m2a(1.0),
stress("stress", id, memory_id),
element_material("element_material", id, memory_id),
set_ID("beam sets", id, memory_id),
stiffness_matrix(NULL),
mass_matrix(NULL),
velocity_damping_matrix(NULL),
jacobian_matrix(NULL),
solver(NULL),
rotation_matrix("rotation_matices", id, memory_id),
increment_flag(false),
integrator(NULL) {
AKANTU_DEBUG_IN();
registerFEEngineObject<MyFEEngineType>("StructuralMechanicsFEEngine", mesh, spatial_dimension);
this->displacement_rotation = NULL;
this->mass = NULL;
this->velocity = NULL;
this->acceleration = NULL;
this->force_momentum = NULL;
this->residual = NULL;
this->blocked_dofs = NULL;
this->increment = NULL;
this->previous_displacement = NULL;
if(spatial_dimension == 2)
nb_degree_of_freedom = 3;
else if (spatial_dimension == 3)
nb_degree_of_freedom = 6;
else {
AKANTU_DEBUG_TO_IMPLEMENT();
}
this->mesh.registerDumper<DumperParaview>("paraview_all", id, true);
this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_structural);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
StructuralMechanicsModel::~StructuralMechanicsModel() {
AKANTU_DEBUG_IN();
delete integrator;
delete solver;
delete stiffness_matrix;
delete jacobian_matrix;
delete mass_matrix;
delete velocity_damping_matrix;
AKANTU_DEBUG_OUT();
}
void StructuralMechanicsModel::setTimeStep(Real time_step) {
this->time_step = time_step;
#if defined(AKANTU_USE_IOHELPER)
this->mesh.getDumper().setTimeStep(time_step);
#endif
}
/* -------------------------------------------------------------------------- */
/* Initialisation */
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::initFull(const ModelOptions & options) {
const StructuralMechanicsModelOptions & stmm_options =
dynamic_cast<const StructuralMechanicsModelOptions &>(options);
method = stmm_options.analysis_method;
initModel();
initArrays();
displacement_rotation->clear();
velocity ->clear();
acceleration ->clear();
force_momentum ->clear();
residual ->clear();
blocked_dofs ->clear();
increment ->clear();
Mesh::type_iterator it = getFEEngine().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural);
Mesh::type_iterator end = getFEEngine().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural);
for (; it != end; ++it) {
computeRotationMatrix(*it);
}
switch(method) {
case _implicit_dynamic:
initImplicit();
break;
case _static:
initSolver();
break;
default:
AKANTU_EXCEPTION("analysis method not recognised by StructuralMechanicsModel");
break;
}
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::initArrays() {
AKANTU_DEBUG_IN();
UInt nb_nodes = getFEEngine().getMesh().getNbNodes();
std::stringstream sstr_disp; sstr_disp << id << ":displacement";
std::stringstream sstr_velo; sstr_velo << id << ":velocity";
std::stringstream sstr_acce; sstr_acce << id << ":acceleration";
std::stringstream sstr_forc; sstr_forc << id << ":force";
std::stringstream sstr_resi; sstr_resi << id << ":residual";
std::stringstream sstr_boun; sstr_boun << id << ":blocked_dofs";
std::stringstream sstr_incr; sstr_incr << id << ":increment";
displacement_rotation = &(alloc<Real>(sstr_disp.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE));
velocity = &(alloc<Real>(sstr_velo.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE));
acceleration = &(alloc<Real>(sstr_acce.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE));
force_momentum = &(alloc<Real>(sstr_forc.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE));
residual = &(alloc<Real>(sstr_resi.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE));
blocked_dofs = &(alloc<bool>(sstr_boun.str(), nb_nodes, nb_degree_of_freedom, false));
increment = &(alloc<Real>(sstr_incr.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE));
Mesh::type_iterator it = getFEEngine().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural);
Mesh::type_iterator end = getFEEngine().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural);
for (; it != end; ++it) {
UInt nb_element = getFEEngine().getMesh().getNbElement(*it);
UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(*it);
element_material.alloc(nb_element, 1, *it, _not_ghost);
set_ID.alloc(nb_element, 1, *it, _not_ghost);
UInt size = getTangentStiffnessVoigtSize(*it);
stress.alloc(nb_element * nb_quadrature_points, size , *it, _not_ghost);
}
dof_synchronizer = new DOFSynchronizer(getFEEngine().getMesh(), nb_degree_of_freedom);
dof_synchronizer->initLocalDOFEquationNumbers();
dof_synchronizer->initGlobalDOFEquationNumbers();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::initModel() {
getFEEngine().initShapeFunctions(_not_ghost);
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::initSolver(__attribute__((unused)) SolverOptions & options) {
AKANTU_DEBUG_IN();
const Mesh & mesh = getFEEngine().getMesh();
#if !defined(AKANTU_USE_MUMPS) // or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake
AKANTU_DEBUG_ERROR("You should at least activate one solver.");
#else
UInt nb_global_node = mesh.getNbGlobalNodes();
std::stringstream sstr; sstr << id << ":jacobian_matrix";
jacobian_matrix = new SparseMatrix(nb_global_node * nb_degree_of_freedom, _symmetric, sstr.str(), memory_id);
jacobian_matrix->buildProfile(mesh, *dof_synchronizer, nb_degree_of_freedom);
std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix";
stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id);
#ifdef AKANTU_USE_MUMPS
std::stringstream sstr_solv; sstr_solv << id << ":solver";
solver = new SolverMumps(*jacobian_matrix, sstr_solv.str());
dof_synchronizer->initScatterGatherCommunicationScheme();
#else
AKANTU_DEBUG_ERROR("You should at least activate one solver.");
#endif //AKANTU_USE_MUMPS
solver->initialize(options);
#endif //AKANTU_HAS_SOLVER
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::initImplicit(bool dynamic, SolverOptions & solver_options) {
AKANTU_DEBUG_IN();
if (!increment) setIncrementFlagOn();
initSolver(solver_options);
if(integrator) delete integrator;
integrator = new TrapezoidalRule2();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
UInt StructuralMechanicsModel::getTangentStiffnessVoigtSize(const ElementType & type) {
UInt size;
#define GET_TANGENT_STIFFNESS_VOIGT_SIZE(type) \
size = getTangentStiffnessVoigtSize<type>();
AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(GET_TANGENT_STIFFNESS_VOIGT_SIZE);
#undef GET_TANGENT_STIFFNESS_VOIGT_SIZE
return size;
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::assembleStiffnessMatrix() {
AKANTU_DEBUG_IN();
stiffness_matrix->clear();
Mesh::type_iterator it = getFEEngine().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural);
Mesh::type_iterator end = getFEEngine().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural);
for (; it != end; ++it) {
ElementType type = *it;
#define ASSEMBLE_STIFFNESS_MATRIX(type) \
assembleStiffnessMatrix<type>();
AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(ASSEMBLE_STIFFNESS_MATRIX);
#undef ASSEMBLE_STIFFNESS_MATRIX
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<>
void StructuralMechanicsModel::computeRotationMatrix<_bernoulli_beam_2>(Array<Real> & rotations){
ElementType type = _bernoulli_beam_2;
Mesh & mesh = getFEEngine().getMesh();
UInt nb_element = mesh.getNbElement(type);
Array<UInt>::iterator< Vector<UInt> > connec_it = mesh.getConnectivity(type).begin(2);
Array<Real>::vector_iterator nodes_it = mesh.getNodes().begin(spatial_dimension);
Array<Real>::matrix_iterator R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom);
for (UInt e = 0; e < nb_element; ++e, ++R_it, ++connec_it) {
Matrix<Real> & R = *R_it;
Vector<UInt> & connec = *connec_it;
Vector<Real> x2;
x2 = nodes_it[connec(1)]; // X2
Vector<Real> x1;
x1 = nodes_it[connec(0)]; // X1
Real le = x1.distance(x2);
Real c = (x2(0) - x1(0)) / le;
Real s = (x2(1) - x1(1)) / le;
/// Definition of the rotation matrix
R(0,0) = c; R(0,1) = s; R(0,2) = 0.;
R(1,0) = -s; R(1,1) = c; R(1,2) = 0.;
R(2,0) = 0.; R(2,1) = 0.; R(2,2) = 1.;
}
}
/* -------------------------------------------------------------------------- */
template<>
void StructuralMechanicsModel::computeRotationMatrix<_bernoulli_beam_3>(Array<Real> & rotations){
ElementType type = _bernoulli_beam_3;
Mesh & mesh = getFEEngine().getMesh();
UInt nb_element = mesh.getNbElement(type);
Array<Real>::vector_iterator n_it = mesh.getNormals(type).begin(spatial_dimension);
Array<UInt>::iterator< Vector<UInt> > connec_it = mesh.getConnectivity(type).begin(2);
Array<Real>::vector_iterator nodes_it = mesh.getNodes().begin(spatial_dimension);
Matrix<Real> Pe (spatial_dimension, spatial_dimension);
Matrix<Real> Pg (spatial_dimension, spatial_dimension);
Matrix<Real> inv_Pg(spatial_dimension, spatial_dimension);
Vector<Real> x_n(spatial_dimension); // x vect n
Array<Real>::matrix_iterator R_it =
rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom);
for (UInt e=0 ; e < nb_element; ++e, ++n_it, ++connec_it, ++R_it) {
Vector<Real> & n = *n_it;
Matrix<Real> & R = *R_it;
Vector<UInt> & connec = *connec_it;
Vector<Real> x;
x = nodes_it[connec(1)]; // X2
Vector<Real> y;
y = nodes_it[connec(0)]; // X1
Real l = x.distance(y);
x -= y; // X2 - X1
x_n.crossProduct(x, n);
Pe.eye();
Pe(0, 0) *= l;
Pe(1, 1) *= -l;
Pg(0,0) = x(0); Pg(0,1) = x_n(0); Pg(0,2) = n(0);
Pg(1,0) = x(1); Pg(1,1) = x_n(1); Pg(1,2) = n(1);
Pg(2,0) = x(2); Pg(2,1) = x_n(2); Pg(2,2) = n(2);
inv_Pg.inverse(Pg);
Pe *= inv_Pg;
for (UInt i = 0; i < spatial_dimension; ++i) {
for (UInt j = 0; j < spatial_dimension; ++j) {
R(i, j) = Pe(i, j);
R(i + spatial_dimension,j + spatial_dimension) = Pe(i, j);
}
}
}
}
/* -------------------------------------------------------------------------- */
template<>
void StructuralMechanicsModel::computeRotationMatrix<_kirchhoff_shell>(Array<Real> & rotations){
ElementType type = _kirchhoff_shell;
Mesh & mesh = getFEEngine().getMesh();
UInt nb_element = mesh.getNbElement(type);
Array<UInt>::iterator< Vector<UInt> > connec_it = mesh.getConnectivity(type).begin(3);
Array<Real>::vector_iterator nodes_it = mesh.getNodes().begin(spatial_dimension);
Matrix<Real> Pe (spatial_dimension, spatial_dimension);
Matrix<Real> Pg (spatial_dimension, spatial_dimension);
Matrix<Real> inv_Pg(spatial_dimension, spatial_dimension);
Array<Real>::matrix_iterator R_it =
rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom);
for (UInt e=0 ; e < nb_element; ++e, ++connec_it, ++R_it) {
Pe.eye();
Matrix<Real> & R = *R_it;
Vector<UInt> & connec = *connec_it;
Vector<Real> x2;
x2 = nodes_it[connec(1)]; // X2
Vector<Real> x1;
x1 = nodes_it[connec(0)]; // X1
Vector<Real> x3;
x3 = nodes_it[connec(2)]; // X3
Vector<Real> Pg_col_1=x2-x1;
Vector<Real> Pg_col_2 = x3-x1;
Vector<Real> Pg_col_3(spatial_dimension);
Pg_col_3.crossProduct(Pg_col_1, Pg_col_2);
for (UInt i = 0; i <spatial_dimension; ++i) {
Pg(i,0) = Pg_col_1(i);
Pg(i,1) = Pg_col_2(i);
Pg(i,2) = Pg_col_3(i);
}
inv_Pg.inverse(Pg);
// Pe *= inv_Pg;
Pe.eye();
for (UInt i = 0; i < spatial_dimension; ++i) {
for (UInt j = 0; j < spatial_dimension; ++j) {
R(i, j) = Pe(i, j);
R(i + spatial_dimension,j + spatial_dimension) = Pe(i, j);
}
}
}
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::computeRotationMatrix(const ElementType & type) {
Mesh & mesh = getFEEngine().getMesh();
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
UInt nb_element = mesh.getNbElement(type);
if(!rotation_matrix.exists(type)) {
rotation_matrix.alloc(nb_element,
nb_degree_of_freedom*nb_nodes_per_element * nb_degree_of_freedom*nb_nodes_per_element,
type);
} else {
rotation_matrix(type).resize(nb_element);
}
rotation_matrix(type).clear();
Array<Real>rotations(nb_element, nb_degree_of_freedom * nb_degree_of_freedom);
rotations.clear();
#define COMPUTE_ROTATION_MATRIX(type) \
computeRotationMatrix<type>(rotations);
AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_ROTATION_MATRIX);
#undef COMPUTE_ROTATION_MATRIX
Array<Real>::matrix_iterator R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom);
Array<Real>::matrix_iterator T_it =
rotation_matrix(type).begin(nb_degree_of_freedom*nb_nodes_per_element,
nb_degree_of_freedom*nb_nodes_per_element);
for (UInt el = 0; el < nb_element; ++el, ++R_it, ++T_it) {
Matrix<Real> & T = *T_it;
Matrix<Real> & R = *R_it;
T.clear();
for (UInt k = 0; k < nb_nodes_per_element; ++k){
for (UInt i = 0; i < nb_degree_of_freedom; ++i)
for (UInt j = 0; j < nb_degree_of_freedom; ++j)
T(k*nb_degree_of_freedom + i, k*nb_degree_of_freedom + j) = R(i, j);
}
}
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::computeStresses() {
AKANTU_DEBUG_IN();
Mesh::type_iterator it = getFEEngine().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural);
Mesh::type_iterator end = getFEEngine().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural);
for (; it != end; ++it) {
ElementType type = *it;
#define COMPUTE_STRESS_ON_QUAD(type) \
computeStressOnQuad<type>();
AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_STRESS_ON_QUAD);
#undef COMPUTE_STRESS_ON_QUAD
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::updateResidual() {
AKANTU_DEBUG_IN();
residual->copy(*force_momentum);
Array<Real> ku(*displacement_rotation, true);
ku *= *stiffness_matrix;
*residual -= ku;
this->updateResidualInternal();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::implicitPred() {
AKANTU_DEBUG_IN();
if(previous_displacement) previous_displacement->copy(*displacement_rotation);
if(method == _implicit_dynamic)
integrator->integrationSchemePred(time_step,
*displacement_rotation,
*velocity,
*acceleration,
*blocked_dofs);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::implicitCorr() {
AKANTU_DEBUG_IN();
Real * incr_val = increment->storage();
bool * boun_val = blocked_dofs->storage();
UInt nb_nodes = displacement_rotation->getSize();
for (UInt j = 0; j < nb_nodes * nb_degree_of_freedom; ++j, ++incr_val, ++boun_val){
*incr_val *= (1. - *boun_val);
}
if(method == _implicit_dynamic) {
integrator->integrationSchemeCorrDispl(time_step,
*displacement_rotation,
*velocity,
*acceleration,
*blocked_dofs,
*increment);
} else {
Real * disp_val = displacement_rotation->storage();
incr_val = increment->storage();
- for (UInt j = 0; j < nb_nodes *nb_degree_of_freedom; ++j, ++disp_val){
+ for (UInt j = 0; j < nb_nodes *nb_degree_of_freedom; ++j, ++disp_val, ++incr_val){
*disp_val += *incr_val;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::updateResidualInternal() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Update the residual");
// f = f_ext - f_int - Ma - Cv = r - Ma - Cv;
if(method != _static) {
// f -= Ma
if(mass_matrix) {
// if full mass_matrix
Array<Real> * Ma = new Array<Real>(*acceleration, true, "Ma");
*Ma *= *mass_matrix;
/// \todo check unit conversion for implicit dynamics
// *Ma /= f_m2a
*residual -= *Ma;
delete Ma;
} else if (mass) {
// else lumped mass
UInt nb_nodes = acceleration->getSize();
UInt nb_degree_of_freedom = acceleration->getNbComponent();
Real * mass_val = mass->storage();
Real * accel_val = acceleration->storage();
Real * res_val = residual->storage();
bool * blocked_dofs_val = blocked_dofs->storage();
for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) {
if(!(*blocked_dofs_val)) {
*res_val -= *accel_val * *mass_val /f_m2a;
}
blocked_dofs_val++;
res_val++;
mass_val++;
accel_val++;
}
} else {
AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix.");
}
// f -= Cv
if(velocity_damping_matrix) {
Array<Real> * Cv = new Array<Real>(*velocity);
*Cv *= *velocity_damping_matrix;
*residual -= *Cv;
delete Cv;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::setIncrementFlagOn() {
AKANTU_DEBUG_IN();
if(!increment) {
UInt nb_nodes = mesh.getNbNodes();
std::stringstream sstr_inc; sstr_inc << id << ":increment";
increment = &(alloc<Real>(sstr_inc.str(), nb_nodes, spatial_dimension, 0.));
}
increment_flag = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::solve() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Solving an implicit step.");
UInt nb_nodes = displacement_rotation->getSize();
/// todo residual = force - Kxr * d_bloq
jacobian_matrix->copyContent(*stiffness_matrix);
jacobian_matrix->applyBoundary(*blocked_dofs);
jacobian_matrix->saveMatrix("Kbound.mtx");
increment->clear();
solver->setRHS(*residual);
solver->factorize();
solver->solve(*increment);
Real * increment_val = increment->storage();
Real * displacement_val = displacement_rotation->storage();
bool * blocked_dofs_val = blocked_dofs->storage();
for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) {
if(!(*blocked_dofs_val)) {
*displacement_val += *increment_val;
}
else {
*increment_val = 0.0;
}
displacement_val++;
blocked_dofs_val++;
increment_val++;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<>
bool StructuralMechanicsModel::testConvergence<_scc_increment>(Real tolerance, Real & error){
AKANTU_DEBUG_IN();
UInt nb_nodes = displacement_rotation->getSize();
UInt nb_degree_of_freedom = displacement_rotation->getNbComponent();
error = 0;
Real norm[2] = {0., 0.};
Real * increment_val = increment->storage();
bool * blocked_dofs_val = blocked_dofs->storage();
Real * displacement_val = displacement_rotation->storage();
for (UInt n = 0; n < nb_nodes; ++n) {
for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
if(!(*blocked_dofs_val)) {
norm[0] += *increment_val * *increment_val;
norm[1] += *displacement_val * *displacement_val;
}
blocked_dofs_val++;
increment_val++;
displacement_val++;
}
}
StaticCommunicator::getStaticCommunicator().allReduce(norm, 2, _so_sum);
norm[0] = sqrt(norm[0]);
norm[1] = sqrt(norm[1]);
AKANTU_DEBUG_ASSERT(!Math::isnan(norm[0]), "Something went wrong in the solve phase");
AKANTU_DEBUG_OUT();
if(norm[1] > Math::getTolerance())
error = norm[0] / norm[1];
else
error = norm[0]; //In case the total displacement is zero!
return (error < tolerance);
}
/* -------------------------------------------------------------------------- */
template<>
bool StructuralMechanicsModel::testConvergence<_scc_residual>(Real tolerance, Real & norm) {
AKANTU_DEBUG_IN();
UInt nb_nodes = residual->getSize();
norm = 0;
Real * residual_val = residual->storage();
bool * blocked_dofs_val = blocked_dofs->storage();
for (UInt n = 0; n < nb_nodes; ++n) {
bool is_local_node = mesh.isLocalOrMasterNode(n);
if(is_local_node) {
- for (UInt d = 0; d < spatial_dimension; ++d) {
+ for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
if(!(*blocked_dofs_val)) {
norm += *residual_val * *residual_val;
}
blocked_dofs_val++;
residual_val++;
}
} else {
blocked_dofs_val += spatial_dimension;
residual_val += spatial_dimension;
}
}
StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum);
norm = sqrt(norm);
AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase");
AKANTU_DEBUG_OUT();
return (norm < tolerance);
}
/* -------------------------------------------------------------------------- */
bool StructuralMechanicsModel::testConvergenceIncrement(Real tolerance) {
Real error;
bool tmp = testConvergenceIncrement(tolerance, error);
AKANTU_DEBUG_INFO("Norm of increment : " << error);
return tmp;
}
/* -------------------------------------------------------------------------- */
bool StructuralMechanicsModel::testConvergenceIncrement(Real tolerance, Real & error) {
AKANTU_DEBUG_IN();
Mesh & mesh= getFEEngine().getMesh();
UInt nb_nodes = displacement_rotation->getSize();
UInt nb_degree_of_freedom = displacement_rotation->getNbComponent();
Real norm = 0;
Real * increment_val = increment->storage();
bool * blocked_dofs_val = blocked_dofs->storage();
for (UInt n = 0; n < nb_nodes; ++n) {
bool is_local_node = mesh.isLocalOrMasterNode(n);
for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
if(!(*blocked_dofs_val) && is_local_node) {
norm += *increment_val * *increment_val;
}
blocked_dofs_val++;
increment_val++;
}
}
StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum);
error = sqrt(norm);
AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase");
AKANTU_DEBUG_OUT();
return (error < tolerance);
}
/* -------------------------------------------------------------------------- */
template<>
void StructuralMechanicsModel::computeTangentModuli<_bernoulli_beam_2>(Array<Real> & tangent_moduli) {
UInt nb_element = getFEEngine().getMesh().getNbElement(_bernoulli_beam_2);
UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(_bernoulli_beam_2);
UInt tangent_size = 2;
Array<Real>::matrix_iterator D_it = tangent_moduli.begin(tangent_size, tangent_size);
Array<UInt> & el_mat = element_material(_bernoulli_beam_2, _not_ghost);
for (UInt e = 0; e < nb_element; ++e) {
UInt mat = el_mat(e);
Real E = materials[mat].E;
Real A = materials[mat].A;
Real I = materials[mat].I;
for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it) {
Matrix<Real> & D = *D_it;
D(0,0) = E * A;
D(1,1) = E * I;
}
}
}
/* -------------------------------------------------------------------------- */
template<>
void StructuralMechanicsModel::computeTangentModuli<_bernoulli_beam_3>(Array<Real> & tangent_moduli) {
UInt nb_element = getFEEngine().getMesh().getNbElement(_bernoulli_beam_3);
UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(_bernoulli_beam_3);
UInt tangent_size = 4;
Array<Real>::matrix_iterator D_it = tangent_moduli.begin(tangent_size, tangent_size);
for (UInt e = 0; e < nb_element; ++e) {
UInt mat = element_material(_bernoulli_beam_3, _not_ghost)(e);
Real E = materials[mat].E;
Real A = materials[mat].A;
Real Iz = materials[mat].Iz;
Real Iy = materials[mat].Iy;
Real GJ = materials[mat].GJ;
for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it) {
Matrix<Real> & D = *D_it;
D(0,0) = E * A;
D(1,1) = E * Iz;
D(2,2) = E * Iy;
D(3,3) = GJ;
}
}
}
/* -------------------------------------------------------------------------- */
template<>
void StructuralMechanicsModel::computeTangentModuli<_kirchhoff_shell>(Array<Real> & tangent_moduli) {
UInt nb_element = getFEEngine().getMesh().getNbElement(_kirchhoff_shell);
UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(_kirchhoff_shell);
UInt tangent_size = 6;
Array<Real>::matrix_iterator D_it = tangent_moduli.begin(tangent_size, tangent_size);
for (UInt e = 0; e < nb_element; ++e) {
UInt mat = element_material(_kirchhoff_shell, _not_ghost)(e);
Real E = materials[mat].E;
Real nu = materials[mat].nu;
Real t = materials[mat].t;
Real HH=E*t/(1-nu*nu);
Real DD=E*t*t*t/(12*(1-nu*nu));
for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it) {
Matrix<Real> & D = *D_it;
D(0,0) = HH;
D(0,1) = HH*nu;
D(1,0) = HH*nu;
D(1,1) = HH;
D(2,2) = HH*(1-nu)/2;
D(3,3) = DD;
D(3,4) = DD*nu;
D(4,3) = DD*nu;
D(4,4) = DD;
D(5,5) = DD*(1-nu)/2;
}
}
}
/* -------------------------------------------------------------------------- */
template<>
void StructuralMechanicsModel::transferBMatrixToSymVoigtBMatrix<_bernoulli_beam_2>(Array<Real> & b, bool local) {
UInt nb_element = getFEEngine().getMesh().getNbElement(_bernoulli_beam_2);
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(_bernoulli_beam_2);
UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(_bernoulli_beam_2);
MyFEEngineType & fem = getFEEngineClass<MyFEEngineType>();
Array<Real>::const_vector_iterator shape_Np = fem.getShapesDerivatives(_bernoulli_beam_2, _not_ghost, 0).begin(nb_nodes_per_element);
Array<Real>::const_vector_iterator shape_Mpp = fem.getShapesDerivatives(_bernoulli_beam_2, _not_ghost, 1).begin(nb_nodes_per_element);
Array<Real>::const_vector_iterator shape_Lpp = fem.getShapesDerivatives(_bernoulli_beam_2, _not_ghost, 2).begin(nb_nodes_per_element);
UInt tangent_size = getTangentStiffnessVoigtSize<_bernoulli_beam_2>();
UInt bt_d_b_size = nb_nodes_per_element * nb_degree_of_freedom;
b.clear();
Array<Real>::matrix_iterator B_it = b.begin(tangent_size, bt_d_b_size);
for (UInt e = 0; e < nb_element; ++e) {
for (UInt q = 0; q < nb_quadrature_points; ++q) {
Matrix<Real> & B = *B_it;
const Vector<Real> & Np = *shape_Np;
const Vector<Real> & Lpp = *shape_Lpp;
const Vector<Real> & Mpp = *shape_Mpp;
B(0,0) = Np(0);
B(0,3) = Np(1);
B(1,1) = Mpp(0);
B(1,2) = Lpp(0);
B(1,4) = Mpp(1);
B(1,5) = Lpp(1);
++B_it;
++shape_Np;
++shape_Mpp;
++shape_Lpp;
}
// ++R_it;
}
}
/* -------------------------------------------------------------------------- */
template<>
void StructuralMechanicsModel::transferBMatrixToSymVoigtBMatrix<_bernoulli_beam_3>(Array<Real> & b, bool local) {
MyFEEngineType & fem = getFEEngineClass<MyFEEngineType>();
UInt nb_element = getFEEngine().getMesh().getNbElement(_bernoulli_beam_3);
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(_bernoulli_beam_3);
UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(_bernoulli_beam_3);
Array<Real>::const_vector_iterator shape_Np = fem.getShapesDerivatives(_bernoulli_beam_3, _not_ghost, 0).begin(nb_nodes_per_element);
Array<Real>::const_vector_iterator shape_Mpp = fem.getShapesDerivatives(_bernoulli_beam_3, _not_ghost, 1).begin(nb_nodes_per_element);
Array<Real>::const_vector_iterator shape_Lpp = fem.getShapesDerivatives(_bernoulli_beam_3, _not_ghost, 2).begin(nb_nodes_per_element);
UInt tangent_size = getTangentStiffnessVoigtSize<_bernoulli_beam_3>();
UInt bt_d_b_size = nb_nodes_per_element * nb_degree_of_freedom;
b.clear();
Array<Real>::matrix_iterator B_it = b.begin(tangent_size, bt_d_b_size);
for (UInt e = 0; e < nb_element; ++e) {
for (UInt q = 0; q < nb_quadrature_points; ++q) {
Matrix<Real> & B = *B_it;
const Vector<Real> & Np = *shape_Np;
const Vector<Real> & Lpp = *shape_Lpp;
const Vector<Real> & Mpp = *shape_Mpp;
B(0,0) = Np(0);
B(0,6) = Np(1);
B(1,1) = Mpp(0);
B(1,5) = Lpp(0);
B(1,7) = Mpp(1);
B(1,11) = Lpp(1);
B(2,2) = Mpp(0);
B(2,4) = -Lpp(0);
B(2,8) = Mpp(1);
B(2,10) = -Lpp(1);
B(3,3) = Np(0);
B(3,9) = Np(1);
++B_it;
++shape_Np;
++shape_Mpp;
++shape_Lpp;
}
}
}
/* -------------------------------------------------------------------------- */
template<>
void StructuralMechanicsModel::transferBMatrixToSymVoigtBMatrix<_kirchhoff_shell>(Array<Real> & b, bool local) {
MyFEEngineType & fem = getFEEngineClass<MyFEEngineType>();
UInt nb_element = getFEEngine().getMesh().getNbElement(_kirchhoff_shell);
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(_kirchhoff_shell);
UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(_kirchhoff_shell);
Array<Real>::const_matrix_iterator shape_Np = fem.getShapesDerivatives(_kirchhoff_shell, _not_ghost, 0).begin(2,nb_nodes_per_element);
Array<Real>::const_matrix_iterator shape_Nx1p = fem.getShapesDerivatives(_kirchhoff_shell, _not_ghost, 1).begin(2,nb_nodes_per_element);
Array<Real>::const_matrix_iterator shape_Nx2p = fem.getShapesDerivatives(_kirchhoff_shell, _not_ghost, 2).begin(2,nb_nodes_per_element);
Array<Real>::const_matrix_iterator shape_Nx3p = fem.getShapesDerivatives(_kirchhoff_shell, _not_ghost, 3).begin(2,nb_nodes_per_element);
Array<Real>::const_matrix_iterator shape_Ny1p = fem.getShapesDerivatives(_kirchhoff_shell, _not_ghost, 4).begin(2,nb_nodes_per_element);
Array<Real>::const_matrix_iterator shape_Ny2p = fem.getShapesDerivatives(_kirchhoff_shell, _not_ghost, 5).begin(2,nb_nodes_per_element);
Array<Real>::const_matrix_iterator shape_Ny3p = fem.getShapesDerivatives(_kirchhoff_shell, _not_ghost, 6).begin(2,nb_nodes_per_element);
UInt tangent_size = getTangentStiffnessVoigtSize<_kirchhoff_shell>();
UInt bt_d_b_size = nb_nodes_per_element * nb_degree_of_freedom;
b.clear();
Array<Real>::matrix_iterator B_it = b.begin(tangent_size, bt_d_b_size);
for (UInt e = 0; e < nb_element; ++e) {
for (UInt q = 0; q < nb_quadrature_points; ++q) {
Matrix<Real> & B = *B_it;
const Matrix<Real> & Np = *shape_Np;
const Matrix<Real> & Nx1p = *shape_Nx1p;
const Matrix<Real> & Nx2p = *shape_Nx2p;
const Matrix<Real> & Nx3p = *shape_Nx3p;
const Matrix<Real> & Ny1p = *shape_Ny1p;
const Matrix<Real> & Ny2p = *shape_Ny2p;
const Matrix<Real> & Ny3p = *shape_Ny3p;
B(0,0) = Np(0,0);
B(0,6) = Np(0,1);
B(0,12) = Np(0,2);
B(1,1) = Np(1,0);
B(1,7) = Np(1,1);
B(1,13) = Np(1,2);
B(2,0) = Np(1,0);
B(2,1) = Np(0,0);
B(2,6) = Np(1,1);
B(2,7) = Np(0,1);
B(2,12) = Np(1,2);
B(2,13) = Np(0,2);
B(3,2) = Nx1p(0,0);
B(3,3) = Nx2p(0,0);
B(3,4) = Nx3p(0,0);
B(3,8) = Nx1p(0,1);
B(3,9) = Nx2p(0,1);
B(3,10) = Nx3p(0,1);
B(3,14) = Nx1p(0,2);
B(3,15) = Nx2p(0,2);
B(3,16) = Nx3p(0,2);
B(4,2) = Ny1p(1,0);
B(4,3) = Ny2p(1,0);
B(4,4) = Ny3p(1,0);
B(4,8) = Ny1p(1,1);
B(4,9) = Ny2p(1,1);
B(4,10) = Ny3p(1,1);
B(4,14) = Ny1p(1,2);
B(4,15) = Ny2p(1,2);
B(4,16) = Ny3p(1,2);
B(5,2) = Nx1p(1,0) + Ny1p(0,0);
B(5,3) = Nx2p(1,0) + Ny2p(0,0);
B(5,4) = Nx3p(1,0) + Ny3p(0,0);
B(5,8) = Nx1p(1,1) + Ny1p(0,1);
B(5,9) = Nx2p(1,1) + Ny2p(0,1);
B(5,10) = Nx3p(1,1) + Ny3p(0,1);
B(5,14) = Nx1p(1,2) + Ny1p(0,2);
B(5,15) = Nx2p(1,2) + Ny2p(0,2);
B(5,16) = Nx3p(1,2) + Ny3p(0,2);
++B_it;
++shape_Np;
++shape_Nx1p;
++shape_Nx2p;
++shape_Nx3p;
++shape_Ny1p;
++shape_Ny2p;
++shape_Ny3p;
}
}
}
/* -------------------------------------------------------------------------- */
dumper::Field * StructuralMechanicsModel
::createNodalFieldBool(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
std::map<std::string,Array<bool>* > uint_nodal_fields;
uint_nodal_fields["blocked_dofs" ] = blocked_dofs;
dumper::Field * field = NULL;
field = mesh.createNodalField(uint_nodal_fields[field_name],group_name);
return field;
}
/* -------------------------------------------------------------------------- */
dumper::Field * StructuralMechanicsModel
::createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
UInt n;
if(spatial_dimension == 2) {
n = 2;
} else n = 3;
dumper::Field * field = NULL;
if (field_name == "displacement"){
field = mesh.createStridedNodalField(displacement_rotation,group_name,n,0,padding_flag);
}
if (field_name == "rotation"){
field = mesh.createStridedNodalField(displacement_rotation,group_name,
nb_degree_of_freedom-n,n,padding_flag);
}
if (field_name == "force"){
field = mesh.createStridedNodalField(force_momentum,group_name,n,0,padding_flag);
}
if (field_name == "momentum"){
field = mesh.createStridedNodalField(force_momentum,group_name,
nb_degree_of_freedom-n,n,padding_flag);
}
if (field_name == "residual"){
field = mesh.createNodalField(residual,group_name,padding_flag);
}
return field;
}
/* -------------------------------------------------------------------------- */
dumper::Field * StructuralMechanicsModel
::createElementalField(const std::string & field_name,
const std::string & group_name,
bool padding_flag,
const ElementKind & kind){
dumper::Field * field = NULL;
if(field_name == "element_index_by_material")
field = mesh.createElementalField<UInt, Vector, dumper::ElementalField >(field_name,group_name,this->spatial_dimension,kind);
return field;
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__