Homec4science
Diffusion akantu (master) Browse

rAKA/fem/#fem.cc#
1d4834b8daf3master

Last Change
Raw File

#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__
c4science · Help