Homec4science
Diffusion akantu (features/eigen) Browse

rAKA/extra_packages/igfem/src/mesh_igfem_spherical_growing_gel_tmpl.hh
033a8225e32afeatures/eigen

Last Change
Raw File

mesh_igfem_spherical_growing_gel_tmpl.hh
View Options

/**
* @file mesh_igfem_spherical_growing_gel.cc
* @author Marco Vocialta <marco.vocialta@epfl.ch>
* @date Mon Sep 21 12:42:11 2015
*
* @brief Implementation of the functions of MeshIgfemSphericalGrowingGel
*
* @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 "mesh_utils.hh"
#include <algorithm>
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template <UInt dim>
MeshIgfemSphericalGrowingGel<dim>::MeshIgfemSphericalGrowingGel(Mesh & mesh)
: mesh(mesh), nb_nodes_fem(mesh.getNbNodes()), nb_enriched_nodes(0),
synchronizer(NULL) {}
/* -------------------------------------------------------------------------- */
template <UInt dim> void MeshIgfemSphericalGrowingGel<dim>::init() {
nb_nodes_fem = mesh.getNbNodes();
for (ghost_type_t::iterator gt = ghost_type_t::begin();
gt != ghost_type_t::end(); ++gt) {
GhostType ghost_type = *gt;
Mesh::type_iterator tit = mesh.firstType(dim, ghost_type);
Mesh::type_iterator tend = mesh.lastType(dim, ghost_type);
for (; tit != tend;
++tit) { // loop to add corresponding IGFEM element types
if (*tit == _triangle_3) {
mesh.addConnectivityType(_igfem_triangle_4, ghost_type);
mesh.addConnectivityType(_igfem_triangle_5, ghost_type);
} else
AKANTU_ERROR("Not ready for mesh type " << *tit);
}
tit = mesh.firstType(dim, ghost_type, _ek_not_defined);
tend = mesh.lastType(dim, ghost_type, _ek_not_defined);
for (; tit != tend; ++tit) {
AKANTU_BOOST_LIST_SWITCH(INSTANTIATOR, ELEMENT_LIST, *tit);
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
void MeshIgfemSphericalGrowingGel<dim>::computeMeshQueryListIntersectionPoint(
const std::list<SK::Sphere_3> & query_list) {
/// store number of currently enriched nodes
this->nb_enriched_nodes = mesh.getNbNodes() - nb_nodes_fem;
UInt nb_old_nodes = mesh.getNbNodes();
for (ghost_type_t::iterator gt = ghost_type_t::begin();
gt != ghost_type_t::end(); ++gt) {
GhostType ghost_type = *gt;
Mesh::type_iterator iit = mesh.firstType(dim, ghost_type, _ek_not_defined);
Mesh::type_iterator iend = mesh.lastType(dim, ghost_type, _ek_not_defined);
for (; iit != iend; ++iit) {
MeshAbstractIntersector<SK::Sphere_3> & intersector =
*intersectors(*iit, ghost_type);
intersector.constructData(ghost_type);
intersector.computeMeshQueryListIntersectionPoint(query_list,
nb_old_nodes);
const Array<Real> intersection_points_current_type =
*(intersector.getIntersectionPoints());
const Array<UInt> & new_node_per_elem = intersector.getNewNodePerElem();
/// Send the new node event
UInt new_points = intersection_points_current_type.getSize();
StaticCommunicator::getStaticCommunicator().allReduce(&new_points, 1,
_so_sum);
if (new_points > 0) {
Array<Real> & nodes = this->mesh.getNodes();
UInt nb_node = nodes.getSize();
Array<Real>::const_vector_iterator intersec_p_it =
intersection_points_current_type.begin(dim);
NewIGFEMNodesEvent new_nodes_event;
for (UInt in = 0; in < intersection_points_current_type.getSize();
++in, ++intersec_p_it) {
nodes.push_back(*intersec_p_it);
new_nodes_event.getList().push_back(nb_node + in);
}
new_nodes_event.setNewNodePerElem(new_node_per_elem);
new_nodes_event.setType(*iit);
new_nodes_event.setGhostType(ghost_type);
this->mesh.sendEvent(new_nodes_event);
}
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
void MeshIgfemSphericalGrowingGel<dim>::removeAdditionalNodes() {
AKANTU_DEBUG_IN();
UInt total_nodes = this->mesh.getNbNodes();
UInt nb_new_enriched_nodes =
total_nodes - this->nb_enriched_nodes - this->nb_nodes_fem;
UInt old_total_nodes = this->nb_nodes_fem + this->nb_enriched_nodes;
UInt total_new_nodes = nb_new_enriched_nodes;
StaticCommunicator::getStaticCommunicator().allReduce(&total_new_nodes, 1,
_so_sum);
if (total_new_nodes == 0)
return;
RemovedNodesEvent remove_nodes(this->mesh);
Array<UInt> & nodes_removed = remove_nodes.getList();
Array<UInt> & new_numbering = remove_nodes.getNewNumbering();
for (UInt nnod = 0; nnod < this->nb_nodes_fem; ++nnod) {
new_numbering(nnod) = nnod;
}
for (UInt nnod = nb_nodes_fem; nnod < old_total_nodes; ++nnod) {
new_numbering(nnod) = UInt(-1);
nodes_removed.push_back(nnod);
}
for (UInt nnod = 0; nnod < nb_new_enriched_nodes; ++nnod) {
new_numbering(nnod + old_total_nodes) = this->nb_nodes_fem + nnod;
}
if (nb_new_enriched_nodes > 0) {
for (UInt gt = _not_ghost; gt <= _ghost; ++gt) {
GhostType ghost_type = (GhostType)gt;
Mesh::type_iterator it =
mesh.firstType(_all_dimensions, ghost_type, _ek_not_defined);
Mesh::type_iterator end =
mesh.lastType(_all_dimensions, ghost_type, _ek_not_defined);
for (; it != end; ++it) {
ElementType type = *it;
Array<UInt> & connectivity_array =
mesh.getConnectivity(type, ghost_type);
UInt nb_nodes_per_element = connectivity_array.getNbComponent();
Array<UInt>::vector_iterator conn_it =
connectivity_array.begin(nb_nodes_per_element);
Array<UInt>::vector_iterator conn_end =
connectivity_array.end(nb_nodes_per_element);
for (; conn_it != conn_end; ++conn_it) {
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
(*conn_it)(n) = new_numbering((*conn_it)(n));
}
}
}
}
}
this->mesh.sendEvent(remove_nodes);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt dim> void MeshIgfemSphericalGrowingGel<dim>::buildIGFEMMesh() {
AKANTU_DEBUG_IN();
NewIGFEMElementsEvent new_elements_event;
UInt total_new_elements = 0;
Array<Element> & new_elements_list = new_elements_event.getList();
Array<Element> & old_elements_list = new_elements_event.getOldElementsList();
RemovedElementsEvent removed_elements_event(this->mesh);
ChangedElementsEvent changed_elements_event(this->mesh);
for (ghost_type_t::iterator gt = ghost_type_t::begin();
gt != ghost_type_t::end(); ++gt) {
GhostType ghost_type = *gt;
mesh.initElementTypeMapArray(removed_elements_event.getNewNumbering(), 1,
dim, ghost_type, false, _ek_not_defined);
mesh.initElementTypeMapArray(changed_elements_event.getNewNumbering(), 1,
dim, ghost_type, false, _ek_not_defined);
/// loop over all types in the mesh for a given ghost type
Mesh::type_iterator iit = mesh.firstType(dim, ghost_type, _ek_not_defined);
Mesh::type_iterator iend = mesh.lastType(dim, ghost_type, _ek_not_defined);
for (; iit != iend; ++iit) {
ElementType type = *iit;
MeshAbstractIntersector<SK::Sphere_3> & intersector =
*intersectors(type, ghost_type);
const Array<UInt> & new_node_per_elem = intersector.getNewNodePerElem();
UInt n_new_el = 0;
Array<UInt> & connectivity = this->mesh.getConnectivity(type, ghost_type);
/// get the connectivities of all types that that may transform
Array<UInt> & connec_igfem_tri_4 =
this->mesh.getConnectivity(_igfem_triangle_4, ghost_type);
Array<UInt> & connec_igfem_tri_5 =
this->mesh.getConnectivity(_igfem_triangle_5, ghost_type);
Array<UInt> & connec_tri_3 =
this->mesh.getConnectivity(_triangle_3, ghost_type);
/// create elements to store the newly generated elements
Element el_tri_3(_triangle_3, 0, ghost_type, _ek_regular);
Element el_igfem_tri_4(_igfem_triangle_4, 0, ghost_type, _ek_igfem);
Element el_igfem_tri5(_igfem_triangle_5, 0, ghost_type, _ek_igfem);
Array<UInt> & new_numbering =
removed_elements_event.getNewNumbering(type, ghost_type);
new_numbering.resize(connectivity.getSize());
/// container for element to be removed
Element old_element(type, 0, ghost_type, Mesh::getKind(type));
for (UInt nel = 0; nel != new_node_per_elem.getSize(); ++nel) {
/// a former IGFEM triangle is transformed into a regular triangle
if ((type != _triangle_3) && (new_node_per_elem(nel, 0) == 0)) {
Vector<UInt> connec_new_elem(3);
connec_new_elem(0) = connectivity(nel, 0);
connec_new_elem(1) = connectivity(nel, 1);
connec_new_elem(2) = connectivity(nel, 2);
/// add the new element in the mesh
UInt nb_triangles_3 = connec_tri_3.getSize();
connec_tri_3.push_back(connec_new_elem);
el_tri_3.element = nb_triangles_3;
new_elements_list.push_back(el_tri_3);
/// the number of the old element in the mesh
old_element.element = nel;
old_elements_list.push_back(old_element);
new_numbering(nel) = UInt(-1);
++n_new_el;
}
/// element is enriched
else if (new_node_per_elem(nel, 0) != 0) {
switch (new_node_per_elem(nel, 0)) {
/// new element is of type igfem_triangle_4
case 1: {
Vector<UInt> connec_new_elem(4);
switch (new_node_per_elem(nel, 2)) {
case 0:
connec_new_elem(0) = connectivity(nel, 2);
connec_new_elem(1) = connectivity(nel, 0);
connec_new_elem(2) = connectivity(nel, 1);
break;
case 1:
connec_new_elem(0) = connectivity(nel, 0);
connec_new_elem(1) = connectivity(nel, 1);
connec_new_elem(2) = connectivity(nel, 2);
break;
case 2:
connec_new_elem(0) = connectivity(nel, 1);
connec_new_elem(1) = connectivity(nel, 2);
connec_new_elem(2) = connectivity(nel, 0);
break;
default:
AKANTU_ERROR("A triangle has only 3 segments not "
<< new_node_per_elem(nel, 2));
break;
}
connec_new_elem(3) = new_node_per_elem(nel, 1);
UInt nb_igfem_triangles_4 = connec_igfem_tri_4.getSize();
connec_igfem_tri_4.push_back(connec_new_elem);
el_igfem_tri_4.element = nb_igfem_triangles_4;
new_elements_list.push_back(el_igfem_tri_4);
break;
}
/// new element is of type igfem_triangle_5
case 2: {
Vector<UInt> connec_new_elem(5);
if ((new_node_per_elem(nel, 2) == 0) &&
(new_node_per_elem(nel, 4) == 1)) {
connec_new_elem(0) = connectivity(nel, 1);
connec_new_elem(1) = connectivity(nel, 2);
connec_new_elem(2) = connectivity(nel, 0);
connec_new_elem(3) = new_node_per_elem(nel, 3);
connec_new_elem(4) = new_node_per_elem(nel, 1);
} else if ((new_node_per_elem(nel, 2) == 0) &&
(new_node_per_elem(nel, 4) == 2)) {
connec_new_elem(0) = connectivity(nel, 0);
connec_new_elem(1) = connectivity(nel, 1);
connec_new_elem(2) = connectivity(nel, 2);
connec_new_elem(3) = new_node_per_elem(nel, 1);
connec_new_elem(4) = new_node_per_elem(nel, 3);
} else if ((new_node_per_elem(nel, 2) == 1) &&
(new_node_per_elem(nel, 4) == 2)) {
connec_new_elem(0) = connectivity(nel, 2);
connec_new_elem(1) = connectivity(nel, 0);
connec_new_elem(2) = connectivity(nel, 1);
connec_new_elem(3) = new_node_per_elem(nel, 3);
connec_new_elem(4) = new_node_per_elem(nel, 1);
} else if ((new_node_per_elem(nel, 2) == 1) &&
(new_node_per_elem(nel, 4) == 0)) {
connec_new_elem(0) = connectivity(nel, 1);
connec_new_elem(1) = connectivity(nel, 2);
connec_new_elem(2) = connectivity(nel, 0);
connec_new_elem(3) = new_node_per_elem(nel, 1);
connec_new_elem(4) = new_node_per_elem(nel, 3);
} else if ((new_node_per_elem(nel, 2) == 2) &&
(new_node_per_elem(nel, 4) == 0)) {
connec_new_elem(0) = connectivity(nel, 0);
connec_new_elem(1) = connectivity(nel, 1);
connec_new_elem(2) = connectivity(nel, 2);
connec_new_elem(3) = new_node_per_elem(nel, 3);
connec_new_elem(4) = new_node_per_elem(nel, 1);
} else if ((new_node_per_elem(nel, 2) == 2) &&
(new_node_per_elem(nel, 4) == 1)) {
connec_new_elem(0) = connectivity(nel, 2);
connec_new_elem(1) = connectivity(nel, 0);
connec_new_elem(2) = connectivity(nel, 1);
connec_new_elem(3) = new_node_per_elem(nel, 1);
connec_new_elem(4) = new_node_per_elem(nel, 3);
} else
AKANTU_ERROR("A triangle has only 3 segments (0 to 2) not "
<< new_node_per_elem(nel, 2) << "and"
<< new_node_per_elem(nel, 4));
UInt nb_igfem_triangles_5 = connec_igfem_tri_5.getSize();
connec_igfem_tri_5.push_back(connec_new_elem);
el_igfem_tri5.element = nb_igfem_triangles_5;
new_elements_list.push_back(el_igfem_tri5);
break;
}
default:
AKANTU_ERROR("IGFEM cannot add " << new_node_per_elem(nel, 0)
<< " nodes to triangles");
break;
}
old_element.element = nel;
old_elements_list.push_back(old_element);
new_numbering(nel) = UInt(-1);
++n_new_el;
}
}
total_new_elements += n_new_el;
}
}
/// update new numbering
for (ghost_type_t::iterator gt = ghost_type_t::begin();
gt != ghost_type_t::end(); ++gt) {
GhostType ghost_type = *gt;
/// loop over all types in the mesh for a given ghost type
Mesh::type_iterator iit = mesh.firstType(dim, ghost_type, _ek_not_defined);
Mesh::type_iterator iend = mesh.lastType(dim, ghost_type, _ek_not_defined);
for (; iit != iend; ++iit) {
ElementType type = *iit;
Array<UInt> & new_numbering =
removed_elements_event.getNewNumbering(type, ghost_type);
UInt el_index = 0;
UInt nb_element = this->mesh.getNbElement(type, ghost_type);
new_numbering.resize(nb_element);
for (UInt e = 0; e < nb_element; ++e) {
if (new_numbering(e) != UInt(-1)) {
new_numbering(e) = el_index;
++el_index;
}
}
changed_elements_event.getNewNumbering(type, ghost_type)
.copy(new_numbering);
}
}
StaticCommunicator::getStaticCommunicator().allReduce(&total_new_elements, 1,
_so_sum);
if (total_new_elements > 0) {
changed_elements_event.getListOld().copy(
new_elements_event.getOldElementsList());
changed_elements_event.getListNew().copy(new_elements_event.getList());
this->mesh.sendEvent(changed_elements_event);
this->mesh.sendEvent(new_elements_event);
Array<Element> & removed_list = removed_elements_event.getList();
removed_list.copy(new_elements_event.getOldElementsList());
this->mesh.sendEvent(removed_elements_event);
}
removeAdditionalNodes();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
void MeshIgfemSphericalGrowingGel<dim>::buildSegmentConnectivityToNodeType() {
Mesh mesh_facets(mesh.initMeshFacets());
MeshUtils::buildSegmentToNodeType(mesh, mesh_facets, synchronizer);
// only the ghost elements are considered
for (UInt g = _not_ghost; g <= _ghost; ++g) {
GhostType ghost_type = (GhostType)g;
Mesh::type_iterator it = mesh_facets.firstType(1, ghost_type);
Mesh::type_iterator end = mesh_facets.lastType(1, ghost_type);
for (; it != end; ++it) {
ElementType type = *it;
const Array<Int> & segment_to_nodetype =
mesh_facets.getData<Int>("segment_to_nodetype", type, ghost_type);
const Array<UInt> & segment_connectivity =
mesh_facets.getConnectivity(type, ghost_type);
// looping over all the segments
Array<UInt>::const_vector_iterator conn_it =
segment_connectivity.begin(segment_connectivity.getNbComponent());
Array<UInt>::const_vector_iterator conn_end =
segment_connectivity.end(segment_connectivity.getNbComponent());
UInt seg_index = 0;
UInt n[2];
for (; conn_it != conn_end; ++conn_it, ++seg_index) {
Int seg_type = segment_to_nodetype(seg_index);
n[0] = (*conn_it)(0);
n[1] = (*conn_it)(1);
if ((mesh.isMasterNode(n[0]) || mesh.isSlaveNode(n[0])) &&
(mesh.isMasterNode(n[1]) || mesh.isSlaveNode(n[1]))) {
std::sort(n, n + 2);
segment_conn_to_node_type[std::make_pair(n[0], n[1])] = seg_type;
}
}
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
void MeshIgfemSphericalGrowingGel<dim>::updateNodeType(
const Array<UInt> & nodes_list, const Array<UInt> & new_node_per_elem,
ElementType type, GhostType ghost_type) {
Array<Int> & nodes_type = mesh.getNodesType();
UInt old_nodes = nodes_type.getSize();
UInt new_nodes = nodes_list.getSize();
// exit this function if the simulation in run in serial
if (old_nodes == 0 || new_nodes == 0)
return;
nodes_type.resize(old_nodes + new_nodes);
Array<bool> checked_node(new_nodes, 1, false);
UInt nb_prim_by_el = 0;
if ((type == _triangle_3) || (type == _igfem_triangle_4) ||
(type == _igfem_triangle_5)) {
nb_prim_by_el = 3;
} else {
AKANTU_ERROR("Not ready for mesh type " << type);
}
// determine the node type for the local, master and slave nodes
const Array<UInt> & connectivity = mesh.getConnectivity(type, ghost_type);
unordered_map<std::pair<UInt, UInt>, Int>::type::iterator seg_type_it;
unordered_map<std::pair<UInt, UInt>, Int>::type::iterator seg_type_end =
segment_conn_to_node_type.end();
for (UInt el = 0; el < new_node_per_elem.getSize(); ++el) {
UInt nb_nodes = new_node_per_elem(el, 0);
for (UInt n = 0; n < nb_nodes; ++n) {
UInt node_index = new_node_per_elem(el, 2 * n + 1);
if (node_index < old_nodes || checked_node(node_index - old_nodes))
continue;
// get the elemental connectivity of the segment associated to the node
UInt segment_index = new_node_per_elem(el, 2 * n + 2);
UInt extreme_nodes[2];
extreme_nodes[0] = segment_index;
extreme_nodes[1] = segment_index + 1;
if (extreme_nodes[1] == nb_prim_by_el)
extreme_nodes[1] = 0;
// get the connectivity of the segment
extreme_nodes[0] = connectivity(el, extreme_nodes[0]);
extreme_nodes[1] = connectivity(el, extreme_nodes[1]);
// if one extreme nodes is pure ghost, then also the new node is pure
// ghost
if (mesh.isPureGhostNode(extreme_nodes[0]) ||
mesh.isPureGhostNode(extreme_nodes[1]))
nodes_type(node_index) = -3;
// if on of the two extreme nodes is local, then also the new node is
// local
else if (mesh.isLocalNode(extreme_nodes[0]) ||
mesh.isLocalNode(extreme_nodes[1]))
nodes_type(node_index) = -1;
// otherwise use the value stored in the map
else {
std::sort(extreme_nodes, extreme_nodes + 2);
seg_type_it = segment_conn_to_node_type.find(
std::make_pair(extreme_nodes[0], extreme_nodes[1]));
AKANTU_DEBUG_ASSERT(seg_type_it != seg_type_end, "Segment not found");
nodes_type(node_index) = seg_type_it->second;
}
checked_node(node_index - old_nodes) = true;
}
}
AKANTU_DEBUG_ASSERT(std::accumulate(checked_node.begin(), checked_node.end(),
0) == Int(checked_node.getSize()),
"Not all new nodes were assigned a node type");
}
__END_AKANTU__
c4science · Help