Homec4science
Diffusion akantu (bugfixes/39-the-python-example-for-phasefield-are-broken-due-to-a-change-to-setintegrationscheme) Browse

rAKA/test/test_model/test_non_local_toolbox/test_pair_computation.cc
0b9711b33a04bugfixes/39-the-python…

Last Change
Raw File

test_pair_computation.cc
View Options

/**
* @file test_pair_computation.cc
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @date Tue Nov 24 10:33:55 2015
*
* @brief test the weight computation with and without grid
*
* @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 "solid_mechanics_model.hh"
#include "test_material_damage.hh"
#include "non_local_manager.hh"
#include "non_local_neighborhood.hh"
#include "dumper_paraview.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
typedef std::vector< std::pair<IntegrationPoint, IntegrationPoint> > PairList;
/* -------------------------------------------------------------------------- */
void computePairs(SolidMechanicsModel & model, PairList * pair_list);
int main(int argc, char *argv[]) {
akantu::initialize("material_remove_damage.dat", argc, argv);
// some configuration variables
const UInt spatial_dimension = 2;
StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator();
Int psize = comm.getNbProc();
Int prank = comm.whoAmI();
// mesh creation and read
Mesh mesh(spatial_dimension);
akantu::MeshPartition * partition = NULL;
if(prank == 0) {
mesh.read("pair_test.msh");
/// partition the mesh
partition = new MeshPartitionScotch(mesh, spatial_dimension);
partition->partitionate(psize);
}
/// model creation
SolidMechanicsModel model(mesh);
model.initParallel(partition);
delete partition;
/// creation of material selector
MeshDataMaterialSelector<std::string> * mat_selector;
mat_selector = new MeshDataMaterialSelector<std::string>("physical_names", model);
model.setMaterialSelector(*mat_selector);
/// model initialization changed to use our material
model.initFull(SolidMechanicsModelOptions(_static, true));
model.registerNewCustomMaterials< TestMaterialDamage<spatial_dimension> >("test_material");
model.initMaterials();
/// dump material index in paraview
model.addDumpField("material_index");
model.dump();
/// compute the pairs by looping over all the quadrature points
PairList pair_list[2];
computePairs(model, pair_list);
NonLocalManager & manager = model.getNonLocalManager();
const PairList * pairs_mat_1 = manager.getNeighborhood("mat_1").getPairLists();
const PairList * pairs_mat_2 = manager.getNeighborhood("mat_2").getPairLists();
/// compare the number of pairs
UInt nb_not_ghost_pairs_grid = pairs_mat_1[0].size() + pairs_mat_2[0].size();
UInt nb_ghost_pairs_grid = pairs_mat_1[1].size() + pairs_mat_2[1].size();
UInt nb_not_ghost_pairs_no_grid = pair_list[0].size();
UInt nb_ghost_pairs_no_grid = pair_list[1].size();
if ((nb_not_ghost_pairs_grid != nb_not_ghost_pairs_no_grid) ||
(nb_ghost_pairs_grid != nb_ghost_pairs_no_grid)) {
std::cout << "The number of pairs is not correct: TEST FAILED!!!" << std::endl;
finalize();
return EXIT_FAILURE;
}
for (UInt i = 0; i < pairs_mat_1[0].size(); ++i) {
std::pair<IntegrationPoint,IntegrationPoint> p = (pairs_mat_1[0])[i];
PairList::const_iterator it = std::find(pair_list[0].begin(), pair_list[0].end(), (pairs_mat_1[0])[i]);
if(it == pair_list[0].end()) {
std::cout << "The pairs are not correct" << std::endl;
finalize();
return EXIT_FAILURE;
}
}
for (UInt i = 0; i < pairs_mat_2[0].size(); ++i) {
std::pair<IntegrationPoint,IntegrationPoint> p = (pairs_mat_2[0])[i];
PairList::const_iterator it = std::find(pair_list[0].begin(), pair_list[0].end(), (pairs_mat_2[0])[i]);
if(it == pair_list[0].end()) {
std::cout << "The pairs are not correct" << std::endl;
finalize();
return EXIT_FAILURE;
}
}
for (UInt i = 0; i < pairs_mat_1[1].size(); ++i) {
std::pair<IntegrationPoint,IntegrationPoint> p = (pairs_mat_1[1])[i];
PairList::const_iterator it = std::find(pair_list[1].begin(), pair_list[1].end(), (pairs_mat_1[1])[i]);
if(it == pair_list[1].end()) {
std::cout << "The pairs are not correct" << std::endl;
finalize();
return EXIT_FAILURE;
}
}
for (UInt i = 0; i < pairs_mat_2[1].size(); ++i) {
std::pair<IntegrationPoint,IntegrationPoint> p = (pairs_mat_2[1])[i];
PairList::const_iterator it = std::find(pair_list[1].begin(), pair_list[1].end(), (pairs_mat_2[1])[i]);
if(it == pair_list[1].end()) {
std::cout << "The pairs are not correct" << std::endl;
finalize();
return EXIT_FAILURE;
}
}
finalize();
return EXIT_SUCCESS;
}
/* -------------------------------------------------------------------------- */
void computePairs(SolidMechanicsModel & model, PairList * pair_list) {
ElementKind kind = _ek_regular;
Mesh & mesh = model.getMesh();
UInt spatial_dimension = model.getSpatialDimension();
/// compute the quadrature points
ElementTypeMapReal quad_coords("quad_coords");
mesh.initElementTypeMapArray(quad_coords, spatial_dimension, spatial_dimension, false, _ek_regular, true);
model.getFEEngine().computeIntegrationPointsCoordinates(quad_coords);
/// loop in a n^2 way over all the quads to generate the pairs
Real neighborhood_radius = 0.5;
Mesh::type_iterator it_1 = mesh.firstType(spatial_dimension, _not_ghost, kind);
Mesh::type_iterator last_type_1 = mesh.lastType(spatial_dimension, _not_ghost, kind);
IntegrationPoint q1;
IntegrationPoint q2;
q1.kind = kind;
q2.kind = kind;
GhostType ghost_type_1 = _not_ghost;
q1.ghost_type = ghost_type_1;
Vector<Real> q1_coords(spatial_dimension);
Vector<Real> q2_coords(spatial_dimension);
for(; it_1 != last_type_1 ; ++it_1) {
ElementType type_1 = *it_1;
q1.type = type_1;
UInt nb_elements_1 = mesh.getNbElement(type_1, ghost_type_1);
UInt nb_quads_1 = model.getFEEngine().getNbIntegrationPoints(type_1);
Array<Real> & quad_coords_1 = quad_coords(q1.type, q1.ghost_type);
Array<Real>::const_vector_iterator coord_it_1 = quad_coords_1.begin(spatial_dimension);
for (UInt e_1 = 0; e_1 < nb_elements_1; ++e_1) {
q1.element = e_1;
UInt mat_index_1 = model.getMaterialByElement(q1.type, q1.ghost_type).begin()[q1.element];
for (UInt q_1 = 0; q_1 < nb_quads_1; ++q_1) {
q1.global_num = nb_quads_1 * e_1 + q_1;
q1.num_point = q_1;
q1_coords = coord_it_1[q1.global_num];
/// loop over all other quads and create pairs for this given quad
for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) {
GhostType ghost_type_2 = *gt;
q2.ghost_type = ghost_type_2;
Mesh::type_iterator it_2 = mesh.firstType(spatial_dimension, ghost_type_2, kind);
Mesh::type_iterator last_type_2 = mesh.lastType(spatial_dimension, ghost_type_2, kind);
for(; it_2 != last_type_2 ; ++it_2) {
ElementType type_2 = *it_2;
q2.type = type_2;
UInt nb_elements_2 = mesh.getNbElement(type_2, ghost_type_2);
UInt nb_quads_2 = model.getFEEngine().getNbIntegrationPoints(type_2);
Array<Real> & quad_coords_2 = quad_coords(q2.type, q2.ghost_type);
Array<Real>::const_vector_iterator coord_it_2 = quad_coords_2.begin(spatial_dimension);
for (UInt e_2 = 0; e_2 < nb_elements_2; ++e_2) {
q2.element = e_2;
UInt mat_index_2 = model.getMaterialByElement(q2.type, q2.ghost_type).begin()[q2.element];
for (UInt q_2 = 0; q_2 < nb_quads_2; ++q_2) {
q2.global_num = nb_quads_2 * e_2 + q_2;
q2.num_point = q_2;
q2_coords = coord_it_2[q2.global_num];
Real distance = q1_coords.distance(q2_coords);
if (mat_index_1 != mat_index_2)
continue;
else if(distance <= neighborhood_radius + Math::getTolerance() &&
(q2.ghost_type == _ghost ||
(q2.ghost_type == _not_ghost && q1.global_num <= q2.global_num))) { // storing only half lists
pair_list[q2.ghost_type].push_back(std::make_pair(q1, q2));
}
}
}
}
}
}
}
}
}
c4science · Help