Homec4science
Diffusion akantu (bugfix/69-segmentation-fault-caused-by-model-applybc-aka-fromtraction-trac-traction) Browse

rAKA/extra_packages/igfem/src/mesh_igfem_spherical_growing_gel_tmpl.hh
607319619283bugfix/69-segmentation…

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_DEBUG_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_DEBUG_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_DEBUG_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_DEBUG_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_DEBUG_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