diff --git a/src/model/solid_mechanics/embedded_interface_model.cc b/src/model/solid_mechanics/embedded_interface_model.cc
index 2f9b4cf3d..cf922d8ab 100644
--- a/src/model/solid_mechanics/embedded_interface_model.cc
+++ b/src/model/solid_mechanics/embedded_interface_model.cc
@@ -1,210 +1,210 @@
/**
* @file embedded_interface_model.cc
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @date creation: Mon Mar 9 2015
* @date last modification: Mon Mar 9 2015
*
* @brief Model of Solid Mechanics with embedded interfaces
*
* @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 "embedded_interface_model.hh"
#include "material_reinforcement.hh"
#ifdef AKANTU_USE_IOHELPER
# include "dumper_paraview.hh"
# include "dumpable_inline_impl.hh"
#endif
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
const EmbeddedInterfaceModelOptions
default_embedded_interface_model_options(_explicit_lumped_mass, false, false);
/* -------------------------------------------------------------------------- */
EmbeddedInterfaceModel::EmbeddedInterfaceModel(Mesh & mesh,
Mesh & primitive_mesh,
UInt spatial_dimension,
const ID & id,
const MemoryID & memory_id) :
SolidMechanicsModel(mesh, spatial_dimension, id, memory_id),
intersector(mesh, primitive_mesh),
interface_mesh(NULL),
primitive_mesh(primitive_mesh),
interface_material_selector(NULL)
{
// This pointer should be deleted by ~SolidMechanicsModel()
MaterialSelector * mat_sel_pointer =
new MeshDataMaterialSelector<std::string>("physical_names", *this);
this->setMaterialSelector(*mat_sel_pointer);
interface_mesh = &(intersector.getInterfaceMesh());
// Create 1D FEEngine on the interface mesh
registerFEEngineObject<MyFEEngineType>("EmbeddedInterfaceFEEngine", *interface_mesh, 1);
}
/* -------------------------------------------------------------------------- */
EmbeddedInterfaceModel::~EmbeddedInterfaceModel() {
delete interface_material_selector;
}
/* -------------------------------------------------------------------------- */
void EmbeddedInterfaceModel::initFull(const ModelOptions & options) {
const EmbeddedInterfaceModelOptions & eim_options =
dynamic_cast<const EmbeddedInterfaceModelOptions &>(options);
// We don't want to initiate materials before shape functions are initialized
SolidMechanicsModelOptions dummy_options(eim_options.analysis_method, true);
// Do no initialize interface_mesh if told so
if (!eim_options.no_init_intersections)
intersector.constructData();
// Initialize interface FEEngine
FEEngine & engine = getFEEngine("EmbeddedInterfaceFEEngine");
engine.initShapeFunctions(_not_ghost);
engine.initShapeFunctions(_ghost);
SolidMechanicsModel::initFull(dummy_options);
// This will call SolidMechanicsModel::initMaterials() last
this->initMaterials();
#if defined(AKANTU_USE_IOHELPER)
this->mesh.registerDumper<DumperParaview>("reinforcement", id);
this->mesh.addDumpMeshToDumper("reinforcement", *interface_mesh,
1, _not_ghost, _ek_regular);
#endif
}
/* -------------------------------------------------------------------------- */
// This function is very similar to SolidMechanicsModel's
void EmbeddedInterfaceModel::initMaterials() {
Element element;
delete interface_material_selector;
interface_material_selector =
new InterfaceMeshDataMaterialSelector<std::string>("physical_names", *this);
for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) {
element.ghost_type = *gt;
Mesh::type_iterator it = interface_mesh->firstType(1, *gt);
Mesh::type_iterator end = interface_mesh->lastType(1, *gt);
for (; it != end ; ++it) {
UInt nb_element = interface_mesh->getNbElement(*it, *gt);
element.type = *it;
Array<UInt> & mat_indexes = material_index.alloc(nb_element, 1, *it, *gt);
for (UInt el = 0 ; el < nb_element ; el++) {
element.element = el;
UInt mat_index = (*interface_material_selector)(element);
AKANTU_DEBUG_ASSERT(mat_index < materials.size(),
"The material selector returned an index that does not exist");
mat_indexes(element.element) = mat_index;
materials.at(mat_index)->addElement(*it, el, *gt);
}
}
}
SolidMechanicsModel::initMaterials();
}
-/**
- * DO NOT REMOVE - This prevents the material reinforcement to register
- * their number of components. Problems arise with AvgHomogenizingFunctor
- * if the material reinforcement gives its number of components for a
- * field. Since AvgHomogenizingFunctor verifies that all the fields
- * have the same number of components, an exception is raised.
- */
-ElementTypeMap<UInt> EmbeddedInterfaceModel::getInternalDataPerElem(const std::string & field_name,
- const ElementKind & kind) {
- if (!(this->isInternal(field_name,kind))) AKANTU_EXCEPTION("unknown internal " << field_name);
-
- for (UInt m = 0; m < materials.size() ; ++m) {
- if (materials[m]->isInternal<Real>(field_name, kind)) {
- Material * mat = NULL;
-
- switch(this->spatial_dimension) {
- case 1:
- mat = dynamic_cast<MaterialReinforcement<1> *>(materials[m]);
- break;
-
- case 2:
- mat = dynamic_cast<MaterialReinforcement<2> *>(materials[m]);
- break;
-
- case 3:
- mat = dynamic_cast<MaterialReinforcement<3> *>(materials[m]);
- break;
- }
-
- if (mat == NULL && field_name != "stress_embedded")
- return materials[m]->getInternalDataPerElem<Real>(field_name,kind);
- else if (mat != NULL && field_name == "stress_embedded")
- return mat->getInternalDataPerElem<Real>(field_name, kind, "EmbeddedInterfaceFEEngine");
- }
- }
-
- return ElementTypeMap<UInt>();
-}
+// /**
+// * DO NOT REMOVE - This prevents the material reinforcement to register
+// * their number of components. Problems arise with AvgHomogenizingFunctor
+// * if the material reinforcement gives its number of components for a
+// * field. Since AvgHomogenizingFunctor verifies that all the fields
+// * have the same number of components, an exception is raised.
+// */
+// ElementTypeMap<UInt> EmbeddedInterfaceModel::getInternalDataPerElem(const std::string & field_name,
+// const ElementKind & kind) {
+// if (!(this->isInternal(field_name,kind))) AKANTU_EXCEPTION("unknown internal " << field_name);
+
+// for (UInt m = 0; m < materials.size() ; ++m) {
+// if (materials[m]->isInternal<Real>(field_name, kind)) {
+// Material * mat = NULL;
+
+// switch(this->spatial_dimension) {
+// case 1:
+// mat = dynamic_cast<MaterialReinforcement<1> *>(materials[m]);
+// break;
+
+// case 2:
+// mat = dynamic_cast<MaterialReinforcement<2> *>(materials[m]);
+// break;
+
+// case 3:
+// mat = dynamic_cast<MaterialReinforcement<3> *>(materials[m]);
+// break;
+// }
+
+// if (mat == NULL && field_name != "stress_embedded")
+// return materials[m]->getInternalDataPerElem<Real>(field_name,kind);
+// else if (mat != NULL && field_name == "stress_embedded")
+// return mat->getInternalDataPerElem<Real>(field_name, kind, "EmbeddedInterfaceFEEngine");
+// }
+// }
+
+// return ElementTypeMap<UInt>();
+// }
/* -------------------------------------------------------------------------- */
void EmbeddedInterfaceModel::addDumpGroupFieldToDumper(const std::string & dumper_name,
const std::string & field_id,
const std::string & group_name,
const ElementKind & element_kind,
bool padding_flag) {
#ifdef AKANTU_USE_IOHELPER
dumper::Field * field = NULL;
// If dumper is reinforcement, create a 1D elemental field
if (dumper_name == "reinforcement")
field = this->createElementalField(field_id, group_name, padding_flag, 1, element_kind);
else {
try {
SolidMechanicsModel::addDumpGroupFieldToDumper(dumper_name, field_id, group_name, element_kind, padding_flag);
} catch (...) {}
}
if (field) {
DumperIOHelper & dumper = mesh.getGroupDumper(dumper_name,group_name);
Model::addDumpGroupFieldToDumper(field_id,field,dumper);
}
#endif
}
__END_AKANTU__
diff --git a/src/model/solid_mechanics/embedded_interface_model.hh b/src/model/solid_mechanics/embedded_interface_model.hh
index 9dc01fc75..26c64b5c7 100644
--- a/src/model/solid_mechanics/embedded_interface_model.hh
+++ b/src/model/solid_mechanics/embedded_interface_model.hh
@@ -1,166 +1,166 @@
/**
* @file embedded_interface_model.hh
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @date creation: Mon Mar 9 2015
* @date last modification: Mon Mar 9 2015
*
* @brief Model of Solid Mechanics with embedded interfaces
*
* @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_EMBEDDED_INTERFACE_MODEL_HH__
#define __AKANTU_EMBEDDED_INTERFACE_MODEL_HH__
#include "aka_common.hh"
#include "solid_mechanics_model.hh"
#include "mesh.hh"
#include "embedded_interface_intersector.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/// Options for the EmbeddedInterfaceModel
struct EmbeddedInterfaceModelOptions : SolidMechanicsModelOptions {
/**
* @brief Constructor for EmbeddedInterfaceModelOptions
* @param analysis_method see SolidMeechanicsModelOptions
* @param no_init_intersections ignores the embedded elements
* @param no_init_materials see SolidMeechanicsModelOptions
*/
EmbeddedInterfaceModelOptions(AnalysisMethod analysis_method = _explicit_lumped_mass,
bool no_init_intersections = false,
bool no_init_materials = false):
SolidMechanicsModelOptions(analysis_method, no_init_materials),
no_init_intersections(no_init_intersections)
{}
/// Ignore the reinforcements
bool no_init_intersections;
};
extern const EmbeddedInterfaceModelOptions default_embedded_interface_model_options;
/**
* @brief Solid mechanics model using the embedded model.
*
* This SolidMechanicsModel subclass implements the embedded model,
* a method used to represent 1D elements in a finite elements model
* (eg. reinforcements in concrete).
*
* In addition to the SolidMechanicsModel properties, this model has
* a mesh of the 1D elements embedded in the model, and an instance of the
* EmbeddedInterfaceIntersector class for the computation of the intersections of the
* 1D elements with the background (bulk) mesh.
*
* @see MaterialReinforcement
*/
class EmbeddedInterfaceModel : public SolidMechanicsModel {
/// Same type as SolidMechanicsModel
typedef FEEngineTemplate<IntegratorGauss, ShapeLagrange, _ek_regular> MyFEEngineType;
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
/**
* @brief Constructor
*
* @param mesh main mesh (concrete)
* @param primitive_mesh mesh of the embedded reinforcement
*/
EmbeddedInterfaceModel(Mesh & mesh,
Mesh & primitive_mesh,
UInt spatial_dimension = _all_dimensions,
const ID & id = "embedded_interface_model",
const MemoryID & memory_id = 0);
/// Destructor
virtual ~EmbeddedInterfaceModel();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// Initialise the model
virtual void initFull(const ModelOptions & options = default_embedded_interface_model_options);
/// Initialise the materials
virtual void initMaterials();
/// Allows filtering of dump fields which need to be dumpes on interface mesh
virtual void addDumpGroupFieldToDumper(const std::string & dumper_name,
const std::string & field_id,
const std::string & group_name,
const ElementKind & element_kind,
bool padding_flag);
- virtual ElementTypeMap<UInt> getInternalDataPerElem(const std::string & field_name,
- const ElementKind & kind);
+ // virtual ElementTypeMap<UInt> getInternalDataPerElem(const std::string & field_name,
+ // const ElementKind & kind);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// Get interface mesh
AKANTU_GET_MACRO(InterfaceMesh, *interface_mesh, Mesh &);
/// Get associated elements
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(InterfaceAssociatedElements,
interface_mesh->getData<Element>("associated_element"),
Element);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// Intersector object to build the interface mesh
EmbeddedInterfaceIntersector intersector;
/// Interface mesh (weak reference)
Mesh * interface_mesh;
/// Mesh used to create the CGAL primitives for intersections
Mesh & primitive_mesh;
/// Material selector for interface
MaterialSelector * interface_material_selector;
};
/// Material selector based on mesh data for interface elements
template<typename T>
class InterfaceMeshDataMaterialSelector : public ElementDataMaterialSelector<T> {
public:
InterfaceMeshDataMaterialSelector(const std::string & name, const EmbeddedInterfaceModel & model, UInt first_index = 1) :
ElementDataMaterialSelector<T>(model.getInterfaceMesh().getData<T>(name), model, first_index)
{}
};
__END_AKANTU__
#endif // __AKANTU_EMBEDDED_INTERFACE_MODEL_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 6466fbd5e..e7f819708 100644
--- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc
+++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc
@@ -1,743 +1,743 @@
/**
* @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,
UInt spatial_dimension,
const Mesh & mesh,
FEEngine & fe_engine,
const ID & id) :
Material(model, dim, mesh, fe_engine, id), // /!\ dim, not spatial_dimension !
model(NULL),
stress_embedded("stress_embedded", *this, 1, fe_engine, this->element_filter),
gradu_embedded("gradu_embedded", *this, 1, fe_engine, this->element_filter),
directing_cosines("directing_cosines", *this, 1, fe_engine, this->element_filter),
pre_stress("pre_stress", *this, 1, fe_engine, this->element_filter),
area(1.0),
shape_derivatives() {
AKANTU_DEBUG_IN();
this->initialize(model);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void MaterialReinforcement<dim>::initialize(SolidMechanicsModel & a_model) {
this->model = dynamic_cast<EmbeddedInterfaceModel *>(&a_model);
AKANTU_DEBUG_ASSERT(this->model != NULL, "MaterialReinforcement needs an EmbeddedInterfaceModel");
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");
// Fool the AvgHomogenizingFunctor
//stress.initialize(dim * dim);
// Reallocate the element filter
this->element_filter.free();
this->model->getInterfaceMesh().initElementTypeMapArray(this->element_filter,
1, 1, false, _ek_regular);
}
/* -------------------------------------------------------------------------- */
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();
stress_embedded.initialize(dim * dim);
gradu_embedded.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();
}
/* -------------------------------------------------------------------------- */
/**
* Background shape derivatives need to be stored per background element
* types but also per embedded element type, which is why they are stored
* in an ElementTypeMap<ElementTypeMapArray<Real> *>. The outer ElementTypeMap
* refers to the embedded types, and the inner refers to the background types.
*/
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(1, *int_ghost_it);
Mesh::type_iterator interface_type_end = interface_mesh.lastType(1, *int_ghost_it);
// Loop over interface types
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);
// Loop over background types
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").getNbIntegrationPoints(int_type);
UInt nb_elements = element_filter(int_type, int_ghost).getSize();
// Alloc the background ElementTypeMapArray
shaped_etma->alloc(nb_elements * nb_quad_points,
dim * nb_points,
back_type);
// Insert the background ElementTypeMapArray in the interface ElementTypeMap
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(1, _not_ghost);
Mesh::type_iterator type_end = interface_mesh.lastType(1, _not_ghost);
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(1, _not_ghost);
Mesh::type_iterator type_end = interface_mesh.lastType(1, _not_ghost);
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").getNbIntegrationPoints(type);
Array<Real> & gradu_vec = gradu_embedded(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 = shape_derivatives(type, ghost_type)->operator()(*back_it, ghost_type);
Array<UInt> * background_filter = new Array<UInt>(nb_element, 1, "background_filter");
filterInterfaceBackgroundElements(*background_filter, *back_it, 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.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 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) {
assembleResidualInterface(type, *type_it, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Computes and assemble the residual. Residual in reinforcement is computed as :
* \f[
* \vec{r} = A_s \int_S{\mathbf{B}^T\mathbf{C}^T \vec{\sigma_s}\,\mathrm{d}s}
* \f]
*/
template<UInt dim>
void MaterialReinforcement<dim>::assembleResidualInterface(const ElementType & interface_type,
const ElementType & background_type,
GhostType ghost_type) {
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, ghost_type);
UInt nodes_per_background_e = Mesh::getNbNodesPerElement(background_type);
UInt nb_quadrature_points = interface_engine.getNbIntegrationPoints(interface_type, ghost_type);
UInt nb_element = elem_filter.getSize();
UInt back_dof = dim * nodes_per_background_e;
Array<Real> & shapesd = shape_derivatives(interface_type, ghost_type)->operator()(background_type, ghost_type);
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.begin(dim, nodes_per_background_e);
Array<Real>::matrix_iterator C_it =
directing_cosines(interface_type, ghost_type).begin(voigt_size, voigt_size);
Array<Real>::matrix_iterator sigma_it =
stress_embedded(interface_type, ghost_type).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;
}
Array<Real> * residual_interface = new Array<Real>(nb_element, back_dof, "residual_interface");
interface_engine.integrate(*integrant,
*residual_interface,
back_dof,
interface_type, ghost_type,
elem_filter);
delete integrant;
Array<UInt> * background_filter = new Array<UInt>(nb_element, 1, "background_filter");
filterInterfaceBackgroundElements(*background_filter,
background_type, interface_type,
ghost_type);
background_engine.assembleArray(*residual_interface, residual,
model->getDOFSynchronizer().getLocalDOFEquationNumbers(),
dim, background_type, ghost_type, *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) {
AKANTU_DEBUG_IN();
filter.resize(0);
filter.clear();
Array<Element> & elements = model->getInterfaceAssociatedElements(interface_type, ghost_type);
Array<UInt> & elem_filter = element_filter(interface_type, 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").getNbIntegrationPoints(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) {
assembleStiffnessMatrixInterface(type, *type_it, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Computes the reinforcement stiffness matrix (Gomes & Awruch, 2001)
* \f[
* \mathbf{K}_e = \sum_{i=1}^R{A_i\int_{S_i}{\mathbf{B}^T
* \mathbf{C}_i^T \mathbf{D}_{s, i} \mathbf{C}_i \mathbf{B}\,\mathrm{d}s}}
* \f]
*/
template<UInt dim>
void MaterialReinforcement<dim>::assembleStiffnessMatrixInterface(const ElementType & interface_type,
const ElementType & background_type,
GhostType ghost_type) {
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, ghost_type);
Array<Real> & grad_u = gradu_embedded(interface_type, ghost_type);
UInt nb_element = elem_filter.getSize();
UInt nodes_per_background_e = Mesh::getNbNodesPerElement(background_type);
UInt nb_quadrature_points = interface_engine.getNbIntegrationPoints(interface_type, ghost_type);
UInt back_dof = dim * nodes_per_background_e;
UInt integrant_size = back_dof;
grad_u.resize(nb_quadrature_points * nb_element);
Array<Real> * tangent_moduli = new Array<Real>(nb_element * nb_quadrature_points,
1, "interface_tangent_moduli");
tangent_moduli->clear();
computeTangentModuli(interface_type, *tangent_moduli, ghost_type);
Array<Real> & shapesd = shape_derivatives(interface_type, ghost_type)->operator()(background_type, ghost_type);
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, ghost_type).begin(voigt_size, voigt_size);
Array<Real>::matrix_iterator B_it =
shapesd.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;
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, ghost_type,
elem_filter);
delete integrant;
Array<UInt> * background_filter = new Array<UInt>(nb_element, 1, "background_filter");
filterInterfaceBackgroundElements(*background_filter,
background_type, interface_type,
ghost_type);
background_engine.assembleMatrix(*K_interface, K, dim, background_type, ghost_type, *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.getNbIntegrationPoints(*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.interpolateOnIntegrationPoints(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);
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();
}
template<UInt dim>
Real MaterialReinforcement<dim>::getEnergy(std::string id) {
AKANTU_DEBUG_IN();
if (id == "potential") {
Real epot = 0.;
computePotentialEnergyByElements();
Mesh::type_iterator
it = element_filter.firstType(spatial_dimension),
end = element_filter.lastType(spatial_dimension);
for (; it != end ; ++it) {
FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine");
epot += interface_engine.integrate(potential_energy(*it, _not_ghost),
*it, _not_ghost,
element_filter(*it, _not_ghost));
epot *= area;
}
return epot;
}
AKANTU_DEBUG_OUT();
return 0;
}
-/* -------------------------------------------------------------------------- */
-template<UInt dim>
-ElementTypeMap<UInt> MaterialReinforcement<dim>::getInternalDataPerElem(const ID & field_name,
- const ElementKind & kind,
- const ID & fe_engine_id) const {
- return Material::getInternalDataPerElem(field_name, kind, "EmbeddedInterfaceFEEngine");
-}
-
-
-/* -------------------------------------------------------------------------- */
-// Author is Guillaume Anciaux, see material.cc
-template<UInt dim>
-void MaterialReinforcement<dim>::flattenInternal(const std::string & field_id,
- ElementTypeMapArray<Real> & internal_flat,
- const GhostType ghost_type,
- ElementKind element_kind) const {
- AKANTU_DEBUG_IN();
-
- if (field_id == "stress_embedded" || field_id == "inelastic_strain") {
- Material::flattenInternalIntern(field_id,
- internal_flat,
- 1,
- ghost_type,
- _ek_not_defined,
- &(this->element_filter),
- &(this->model->getInterfaceMesh()));
- }
-
- AKANTU_DEBUG_OUT();
-}
+// /* -------------------------------------------------------------------------- */
+// template<UInt dim>
+// ElementTypeMap<UInt> MaterialReinforcement<dim>::getInternalDataPerElem(const ID & field_name,
+// const ElementKind & kind,
+// const ID & fe_engine_id) const {
+// return Material::getInternalDataPerElem(field_name, kind, "EmbeddedInterfaceFEEngine");
+// }
+
+
+// /* -------------------------------------------------------------------------- */
+// // Author is Guillaume Anciaux, see material.cc
+// template<UInt dim>
+// void MaterialReinforcement<dim>::flattenInternal(const std::string & field_id,
+// ElementTypeMapArray<Real> & internal_flat,
+// const GhostType ghost_type,
+// ElementKind element_kind) const {
+// AKANTU_DEBUG_IN();
+
+// if (field_id == "stress_embedded" || field_id == "inelastic_strain") {
+// Material::flattenInternalIntern(field_id,
+// internal_flat,
+// 1,
+// ghost_type,
+// _ek_not_defined,
+// &(this->element_filter),
+// &(this->model->getInterfaceMesh()));
+// }
+
+// AKANTU_DEBUG_OUT();
+// }
/* -------------------------------------------------------------------------- */
INSTANTIATE_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 e6b5ffeb6..1969cd56f 100644
--- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh
+++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh
@@ -1,197 +1,197 @@
/**
* @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__
/**
* @brief Material used to represent embedded reinforcements
*
* This class is used for computing the reinforcement stiffness matrix
* along with the reinforcement residual. Room is made for constitutive law,
* but actual use of contitutive laws is made in MaterialReinforcementTemplate.
*
* Be careful with the dimensions in this class :
* - this->spatial_dimension is always 1
* - the template parameter dim is the dimension of the problem
*/
template<UInt dim>
class MaterialReinforcement : virtual public Material {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
/// Constructor
MaterialReinforcement(SolidMechanicsModel & model,
UInt spatial_dimension,
const Mesh & mesh,
FEEngine & fe_engine,
const ID & id = "");
/// Destructor
virtual ~MaterialReinforcement();
protected:
void initialize(SolidMechanicsModel & a_model);
/* ------------------------------------------------------------------------ */
/* 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;
/// Compute energy
virtual Real getEnergy(std::string id);
- virtual ElementTypeMap<UInt> getInternalDataPerElem(const ID & field_name,
- const ElementKind & kind,
- const ID & fe_engine_id) const;
-
- /// Reimplementation of Material's function to accomodate for interface mesh
- virtual void flattenInternal(const std::string & field_id,
- ElementTypeMapArray<Real> & internal_flat,
- const GhostType ghost_type = _not_ghost,
- ElementKind element_kind = _ek_not_defined) const;
+ // virtual ElementTypeMap<UInt> getInternalDataPerElem(const ID & field_name,
+ // const ElementKind & kind,
+ // const ID & fe_engine_id) const;
+
+ // /// Reimplementation of Material's function to accomodate for interface mesh
+ // virtual void flattenInternal(const std::string & field_id,
+ // ElementTypeMapArray<Real> & internal_flat,
+ // const GhostType ghost_type = _not_ghost,
+ // ElementKind element_kind = _ek_not_defined) const;
/* ------------------------------------------------------------------------ */
/* 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 assembleStiffnessMatrixInterface(const ElementType & interface_type,
const ElementType & background_type,
GhostType ghost_type);
/// Compute the background shape derivatives for a type
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);
/// 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 assembleResidualInterface(const ElementType & interface_type,
const ElementType & background_type,
GhostType ghost_type);
// 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;
/// Stress in the reinforcement
InternalField<Real> stress_embedded;
/// Gradu of concrete on reinforcement
InternalField<Real> gradu_embedded;
/// C matrix on quad
InternalField<Real> directing_cosines;
/// Prestress on quad
InternalField<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__