material_reinforcement.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Fri, Aug 9, 13:12

material_reinforcement.cc
View Options

	/**
	* @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");


	// 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").getNbQuadraturePoints(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").getNbQuadraturePoints(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.getNbQuadraturePoints(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").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) {
	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.getNbQuadraturePoints(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^tC^tDCB");
	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.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);

	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();

	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__

material_reinforcement.ccNo OneTemporaryActions

File Metadata

material_reinforcement.ccView Options

Event Timeline

material_reinforcement.cc
No OneTemporary
Actions

material_reinforcement.cc
View Options