#fem.cc#
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Mon, Jul 1, 09:18

#fem.cc#
View Options

	/**
	* @file fem.cc
	* @author Nicolas Richart <nicolas.richart@epfl.ch>
	* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
	* @date Fri Jul 16 11:03:02 2010
	*
	* @brief Implementation of the FEM class
	*
	* @section LICENSE
	*
	* Copyright (©) 2010-2011 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 "fem.hh"
	#include "mesh.hh"
	#include "element_class.hh"
	#include "aka_math.hh"


	/* -------------------------------------------------------------------------- */

	__BEGIN_AKANTU__


	/* -------------------------------------------------------------------------- */
	FEM::FEM(Mesh & mesh, UInt element_dimension, FEMID id, MemoryID memory_id) :
	Memory(memory_id), id(id) {
	AKANTU_DEBUG_IN();
	this->element_dimension = (element_dimension != 0) ?
	element_dimension : mesh.getSpatialDimension();

	init();

	this->mesh = &mesh;
	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void FEM::init() {
	for(UInt t = _not_defined; t < _max_element_type; ++t) {
	this->normals_on_quad_points[t] = NULL;
	}
	}

	/* -------------------------------------------------------------------------- */
	FEM::~FEM() {
	AKANTU_DEBUG_IN();

	mesh = NULL;

	AKANTU_DEBUG_OUT();
	}




	/* -------------------------------------------------------------------------- */
	<<<<<<< .mine
	=======
	void FEM::computeNormalsOnQuadPoints(GhostType ghost_type) {
	AKANTU_DEBUG_IN();

	Real * coord = mesh->getNodes().values;
	UInt spatial_dimension = mesh->getSpatialDimension();

	//allocate the normal arrays
	if (ghost_type == _not_ghost)
	mesh->initByElementTypeRealVector(normals_on_quad_points, spatial_dimension, element_dimension,
	id,"normals_onquad",ghost_type);
	else{
	AKANTU_DEBUG_ERROR("to be implemented");
	}

	//loop over the type to build the normals
	const Mesh::ConnectivityTypeList & type_list = mesh->getConnectivityTypeList();
	Mesh::ConnectivityTypeList::const_iterator it;

	for(it = type_list.begin();
	it != type_list.end();
	++it) {

	ElementType type = *it;

	UInt element_type_spatial_dimension = Mesh::getSpatialDimension(type);
	UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
	UInt nb_quad_points = FEM::getNbQuadraturePoints(type);

	if(element_type_spatial_dimension != element_dimension) continue;

	UInt * elem_val;
	UInt nb_element;

	Real * normals_on_quad_val = NULL;

	if(ghost_type == _not_ghost) {
	elem_val = mesh->getConnectivity(type).values;
	nb_element = mesh->getConnectivity(type).getSize();

	normals_on_quad_points[type]->resize(nb_element * nb_quad_points);
	normals_on_quad_val = normals_on_quad_points[type]->values;
	} else {
	elem_val = mesh->getGhostConnectivity(type).values;
	nb_element = mesh->getGhostConnectivity(type).getSize();
	}


	/* ---------------------------------------------------------------------- */
	#define COMPUTE_NORMALS_ON_QUAD(type) \
	do { \
	Real local_coord[spatial_dimension * nb_nodes_per_element]; \
	for (UInt elem = 0; elem < nb_element; ++elem) { \
	mesh->extractNodalCoordinatesFromElement(local_coord, \
	coord,elem_val+elem*nb_nodes_per_element,nb_nodes_per_element); \
	ElementClass<type>::computeNormalsOnQuadPoint(local_coord, \
	spatial_dimension, \
	normals_on_quad_val); \
	normals_on_quad_val += spatial_dimension*nb_quad_points; \
	} \
	} while(0)
	/* ---------------------------------------------------------------------- */

	AKANTU_BOOST_ELEMENT_SWITCH(COMPUTE_NORMALS_ON_QUAD)
	#undef COMPUTE_NORMALS_ON_QUAD

	}
	AKANTU_DEBUG_OUT();
	}


	/* -------------------------------------------------------------------------- */
	void FEM::interpolateOnQuadraturePoints(const Vector<Real> &in_u,
	Vector<Real> &out_uq,
	UInt nb_degre_of_freedom,
	const ElementType & type,
	GhostType ghost_type,
	const Vector<UInt> * filter_elements) const {
	AKANTU_DEBUG_IN();

	Vector<Real> * shapes_loc;
	UInt nb_element;
	UInt * conn_val;

	if(ghost_type == _not_ghost) {
	shapes_loc = shapes[type];
	nb_element = mesh->getNbElement(type);
	conn_val = mesh->getConnectivity(type).values;
	} else {
	shapes_loc = ghost_shapes[type];
	nb_element = mesh->getNbGhostElement(type);
	conn_val = mesh->getGhostConnectivity(type).values;
	}

	AKANTU_DEBUG_ASSERT(shapes_loc != NULL,
	"No shapes for the type " << type);

	UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
	UInt nb_quadrature_points = FEM::getNbQuadraturePoints(type);
	UInt size_of_shapes = FEM::getShapeSize(type);

	AKANTU_DEBUG_ASSERT(in_u.getSize() == mesh->getNbNodes(),
	"The vector in_u(" << in_u.getID()
	<< ") has not the good size.");
	AKANTU_DEBUG_ASSERT(in_u.getNbComponent() == nb_degre_of_freedom,
	"The vector in_u(" << in_u.getID()
	<< ") has not the good number of component.");

	AKANTU_DEBUG_ASSERT(out_uq.getNbComponent() == nb_degre_of_freedom ,
	"The vector out_uq(" << out_uq.getID()
	<< ") has not the good number of component.");

	UInt * filter_elem_val = NULL;
	if(filter_elements != NULL) {
	nb_element = filter_elements->getSize();
	filter_elem_val = filter_elements->values;
	}

	out_uq.resize(nb_element * nb_quadrature_points);

	Real * shape_val = shapes_loc->values;
	Real * u_val = in_u.values;
	Real * uq_val = out_uq.values;

	UInt offset_uq = out_uq.getNbComponent()*nb_quadrature_points;

	Real * shape = shape_val;
	Real * u = static_cast<Real >(calloc(nb_nodes_per_element nb_degre_of_freedom,
	sizeof(Real)));

	for (UInt el = 0; el < nb_element; ++el) {
	UInt el_offset = el * nb_nodes_per_element;
	if(filter_elements != NULL) {
	shape = shape_val + filter_elem_val[el] * size_of_shapes*nb_quadrature_points;
	el_offset = filter_elem_val[el] * nb_nodes_per_element;
	}

	for (UInt n = 0; n < nb_nodes_per_element; ++n) {
	memcpy(u + n * nb_degre_of_freedom,
	u_val + conn_val[el_offset + n] * nb_degre_of_freedom,
	nb_degre_of_freedom * sizeof(Real));
	}

	/// Uq = Shape * U : matrix product
	Math::matrix_matrix(nb_quadrature_points, nb_degre_of_freedom, nb_nodes_per_element,
	shape, u, uq_val);

	uq_val += offset_uq;
	if(filter_elements == NULL) {
	shape += size_of_shapes*nb_quadrature_points;
	}
	}

	free(u);

	#undef INIT_VARIABLES
	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void FEM::gradientOnQuadraturePoints(const Vector<Real> &in_u,
	Vector<Real> &out_nablauq,
	UInt nb_degre_of_freedom,
	const ElementType & type,
	GhostType ghost_type,
	const Vector<UInt> * filter_elements) const {
	AKANTU_DEBUG_IN();

	Vector<Real> * shapesd_loc;
	UInt nb_element;
	UInt * conn_val;

	if(ghost_type == _not_ghost) {
	shapesd_loc = shapes_derivatives[type];
	nb_element = mesh->getNbElement(type);
	conn_val = mesh->getConnectivity(type).values;
	} else {
	shapesd_loc = ghost_shapes_derivatives[type];
	nb_element = mesh->getNbGhostElement(type);
	conn_val = mesh->getGhostConnectivity(type).values;
	}

	AKANTU_DEBUG_ASSERT(shapesd_loc != NULL,
	"No shapes for the type " << type);

	UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
	UInt size_of_shapes_derivatives = FEM::getShapeDerivativesSize(type);
	UInt nb_quadrature_points = FEM::getNbQuadraturePoints(type);

	UInt * filter_elem_val = NULL;
	if(filter_elements != NULL) {
	nb_element = filter_elements->getSize();
	filter_elem_val = filter_elements->values;
	}

	AKANTU_DEBUG_ASSERT(in_u.getSize() == mesh->getNbNodes(),
	"The vector in_u(" << in_u.getID()
	<< ") has not the good size.");
	AKANTU_DEBUG_ASSERT(in_u.getNbComponent() == nb_degre_of_freedom ,
	"The vector in_u(" << in_u.getID()
	<< ") has not the good number of component.");

	AKANTU_DEBUG_ASSERT(out_nablauq.getNbComponent()
	== nb_degre_of_freedom * element_dimension,
	"The vector out_nablauq(" << out_nablauq.getID()
	<< ") has not the good number of component.");

	out_nablauq.resize(nb_element * nb_quadrature_points);

	Real * shaped_val = shapesd_loc->values;
	Real * u_val = in_u.values;
	Real * nablauq_val = out_nablauq.values;

	UInt offset_nablauq = nb_degre_of_freedom * element_dimension;
	UInt offset_shaped = nb_nodes_per_element * element_dimension;

	Real * shaped = shaped_val;
	Real * u = static_cast<Real >(calloc(nb_nodes_per_element nb_degre_of_freedom,
	sizeof(Real)));

	for (UInt el = 0; el < nb_element; ++el) {
	UInt el_offset = el * nb_nodes_per_element;
	if(filter_elements != NULL) {
	shaped = shaped_val + filter_elem_val[el] * size_of_shapes_derivatives*nb_quadrature_points;
	el_offset = filter_elem_val[el] * nb_nodes_per_element;
	}

	for (UInt n = 0; n < nb_nodes_per_element; ++n) {
	memcpy(u + n * nb_degre_of_freedom,
	u_val + conn_val[el_offset + n] * nb_degre_of_freedom,
	nb_degre_of_freedom * sizeof(Real));
	}

	for (UInt q = 0; q < nb_quadrature_points; ++q) {
	/// \nabla(U) = U^t * dphi/dx
	Math::matrixt_matrix(nb_degre_of_freedom, element_dimension, nb_nodes_per_element,
	u,
	shaped,
	nablauq_val);

	nablauq_val += offset_nablauq;
	shaped += offset_shaped;
	}
	}

	free(u);

	#undef INIT_VARIABLES
	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void FEM::integrate(const Vector<Real> & in_f,
	Vector<Real> &intf,
	UInt nb_degre_of_freedom,
	const ElementType & type,
	GhostType ghost_type,
	const Vector<UInt> * filter_elements) const {
	AKANTU_DEBUG_IN();

	Vector<Real> * jac_loc;
	UInt nb_element;

	if(ghost_type == _not_ghost) {
	jac_loc = jacobians[type];
	nb_element = mesh->getNbElement(type);
	} else {
	jac_loc = ghost_jacobians[type];
	nb_element = mesh->getNbGhostElement(type);
	}

	UInt nb_quadrature_points = FEM::getNbQuadraturePoints(type);

	UInt * filter_elem_val = NULL;
	if(filter_elements != NULL) {
	nb_element = filter_elements->getSize();
	filter_elem_val = filter_elements->values;
	}

	AKANTU_DEBUG_ASSERT(in_f.getSize() == nb_element * nb_quadrature_points,
	"The vector in_f(" << in_f.getID() << " size " << in_f.getSize()
	<< ") has not the good size (" << nb_element << ").");
	AKANTU_DEBUG_ASSERT(in_f.getNbComponent() == nb_degre_of_freedom ,
	"The vector in_f(" << in_f.getID()
	<< ") has not the good number of component.");
	AKANTU_DEBUG_ASSERT(intf.getNbComponent() == nb_degre_of_freedom,
	"The vector intf(" << intf.getID()
	<< ") has not the good number of component.");


	intf.resize(nb_element);

	Real * in_f_val = in_f.values;
	Real * intf_val = intf.values;
	Real * jac_val = jac_loc->values;

	UInt offset_in_f = in_f.getNbComponent()*nb_quadrature_points;
	UInt offset_intf = intf.getNbComponent();

	Real * jac = jac_val;

	for (UInt el = 0; el < nb_element; ++el) {
	if(filter_elements != NULL) {
	jac = jac_val + filter_elem_val[el] * nb_quadrature_points;
	}

	integrate(in_f_val, jac, intf_val, nb_degre_of_freedom, nb_quadrature_points);

	in_f_val += offset_in_f;
	intf_val += offset_intf;
	if(filter_elements == NULL) {
	jac += nb_quadrature_points;
	}
	}

	AKANTU_DEBUG_OUT();
	}


	/* -------------------------------------------------------------------------- */
	Real FEM::integrate(const Vector<Real> & in_f,
	const ElementType & type,
	GhostType ghost_type,
	const Vector<UInt> * filter_elements) const {
	AKANTU_DEBUG_IN();
	Vector<Real> * jac_loc;
	UInt nb_element;

	if(ghost_type == _not_ghost) {
	jac_loc = jacobians[type];
	nb_element = mesh->getNbElement(type);
	} else {
	jac_loc = ghost_jacobians[type];
	nb_element = mesh->getNbGhostElement(type);
	}

	UInt nb_quadrature_points = FEM::getNbQuadraturePoints(type);

	UInt * filter_elem_val = NULL;
	if(filter_elements != NULL) {
	nb_element = filter_elements->getSize();
	filter_elem_val = filter_elements->values;
	}

	AKANTU_DEBUG_ASSERT(in_f.getSize() == nb_element * nb_quadrature_points,
	"The vector in_f(" << in_f.getID()
	<< ") has not the good size.");
	AKANTU_DEBUG_ASSERT(in_f.getNbComponent() == 1,
	"The vector in_f(" << in_f.getID()
	<< ") has not the good number of component.");

	Real intf = 0.;
	Real * in_f_val = in_f.values;
	Real * jac_val = jac_loc->values;
	UInt offset_in_f = in_f.getNbComponent() * nb_quadrature_points;
	Real * jac = jac_val;

	for (UInt el = 0; el < nb_element; ++el) {
	if(filter_elements != NULL) {
	jac = jac_val + filter_elem_val[el] * nb_quadrature_points;
	}

	Real el_intf = 0;
	integrate(in_f_val, jac, &el_intf, 1, nb_quadrature_points);
	intf += el_intf;

	in_f_val += offset_in_f;
	if(filter_elements == NULL) {
	jac += nb_quadrature_points;
	}
	}

	AKANTU_DEBUG_OUT();
	return intf;
	}

	/* -------------------------------------------------------------------------- */
	>>>>>>> .r2532
	void FEM::assembleVector(const Vector<Real> & elementary_vect,
	Vector<Real> & nodal_values,
	UInt nb_degre_of_freedom,
	const ElementType & type,
	GhostType ghost_type,
	const Vector<UInt> * filter_elements,
	Real scale_factor) const {
	AKANTU_DEBUG_IN();

	UInt nb_element;
	UInt * conn_val;

	if(ghost_type == _not_ghost) {
	nb_element = mesh->getNbElement(type);
	conn_val = mesh->getConnectivity(type).values;
	} else {
	nb_element = mesh->getNbGhostElement(type);
	conn_val = mesh->getGhostConnectivity(type).values;
	}

	UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
	UInt nb_nodes = mesh->getNbNodes();

	UInt * filter_elem_val = NULL;
	if(filter_elements != NULL) {
	nb_element = filter_elements->getSize();
	filter_elem_val = filter_elements->values;
	}

	AKANTU_DEBUG_ASSERT(elementary_vect.getSize() == nb_element,
	"The vector elementary_vect(" << elementary_vect.getID()
	<< ") has not the good size.");

	AKANTU_DEBUG_ASSERT(elementary_vect.getNbComponent()
	== nb_degre_of_freedom*nb_nodes_per_element,
	"The vector elementary_vect(" << elementary_vect.getID()
	<< ") has not the good number of component.");

	AKANTU_DEBUG_ASSERT(nodal_values.getNbComponent() == nb_degre_of_freedom,
	"The vector nodal_values(" << nodal_values.getID()
	<< ") has not the good number of component.");

	nodal_values.resize(nb_nodes);

	Real * elementary_vect_val = elementary_vect.values;
	Real * nodal_values_val = nodal_values.values;

	for (UInt el = 0; el < nb_element; ++el) {
	UInt el_offset = el * nb_nodes_per_element;
	if(filter_elements != NULL) {
	el_offset = filter_elem_val[el] * nb_nodes_per_element;
	}
	for (UInt n = 0; n < nb_nodes_per_element; ++n) {
	UInt node = conn_val[el_offset + n];
	UInt offset_node = node * nb_degre_of_freedom;
	for (UInt d = 0; d < nb_degre_of_freedom; ++d) {
	nodal_values_val[offset_node + d] += scale_factor * elementary_vect_val[d];
	}
	elementary_vect_val += nb_degre_of_freedom;
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void FEM::printself(std::ostream & stream, int indent) const {
	std::string space;
	for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);

	stream << space << "FEM [" << std::endl;
	stream << space << " + id : " << id << std::endl;
	stream << space << " + element dimension : " << element_dimension << std::endl;

	stream << space << " + mesh [" << std::endl;
	mesh->printself(stream, indent + 2);
	stream << space << AKANTU_INDENT << "]" << std::endl;

	stream << space << " + connectivity type information [" << std::endl;
	const Mesh::ConnectivityTypeList & type_list = mesh->getConnectivityTypeList();
	Mesh::ConnectivityTypeList::const_iterator it;
	for(it = type_list.begin(); it != type_list.end(); ++it) {
	if (mesh->getSpatialDimension(*it) != element_dimension) continue;
	stream << space << AKANTU_INDENT << AKANTU_INDENT << " + " << *it <<" [" << std::endl;
	// if(shapes[*it]) {
	// shapes [*it]->printself(stream, indent + 3);
	// shapes_derivatives[*it]->printself(stream, indent + 3);
	// jacobians [*it]->printself(stream, indent + 3);
	// }
	stream << space << AKANTU_INDENT << AKANTU_INDENT << "]" << std::endl;
	}
	stream << space << AKANTU_INDENT << "]" << std::endl;

	stream << space << "]" << std::endl;
	}


	__END_AKANTU__

#fem.cc#No OneTemporaryActions

File Metadata

#fem.cc#View Options

Event Timeline

#fem.cc#
No OneTemporary
Actions

#fem.cc#
View Options