diff --git a/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc b/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc index 0c480944e..dd1dd08b8 100644 --- a/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc +++ b/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc @@ -1,551 +1,548 @@ /** * @file material_non_local_inline_impl.cc * @author Nicolas Richart * @date Thu Aug 25 11:59:39 2011 * * @brief * * @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 . * */ /* -------------------------------------------------------------------------- */ __END_AKANTU__ /* -------------------------------------------------------------------------- */ #include "aka_types.hh" #include "grid_synchronizer.hh" #include "synchronizer_registry.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialNonLocal::MaterialNonLocal(Model & model, const ID & id) : Material(model, id), weight_func(NULL), cell_list(NULL), update_weigths(0), compute_stress_calls(0) { AKANTU_DEBUG_IN(); is_non_local = true; this->weight_func = new WeightFunction(*this); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template MaterialNonLocal::~MaterialNonLocal() { AKANTU_DEBUG_IN(); delete cell_list; delete weight_func; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::initMaterial() { AKANTU_DEBUG_IN(); // Material::initMaterial(); ByElementTypeReal quadrature_points_coordinates("quadrature_points_coordinates_tmp_nl", id); computeQuadraturePointsCoordinates(quadrature_points_coordinates); weight_func->setRadius(radius); weight_func->init(); createCellList(quadrature_points_coordinates); updatePairList(quadrature_points_coordinates); computeWeights(quadrature_points_coordinates); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::updateResidual(Vector & displacement, GhostType ghost_type) { AKANTU_DEBUG_IN(); // Update the weights for the non local variable averaging if(ghost_type == _not_ghost && this->update_weigths && (this->compute_stress_calls % this->update_weigths == 0)) { ByElementTypeReal quadrature_points_coordinates("quadrature_points_coordinates", id); computeQuadraturePointsCoordinates(quadrature_points_coordinates); computeWeights(quadrature_points_coordinates); } if(ghost_type == _not_ghost) ++this->compute_stress_calls; computeStress(displacement, ghost_type); computeNonLocalStress(ghost_type); assembleResidual(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void MaterialNonLocal::weightedAvergageOnNeighbours(const ByElementTypeVector & to_accumulate, ByElementTypeVector & accumulated, UInt nb_degree_of_freedom) const { AKANTU_DEBUG_IN(); std::set< std::pair >::const_iterator first_pair_types = existing_pairs.begin(); std::set< std::pair >::const_iterator last_pair_types = existing_pairs.end(); GhostType ghost_type1, ghost_type2; ghost_type1 = ghost_type2 = _not_ghost; for (; first_pair_types != last_pair_types; ++first_pair_types) { const Vector & pairs = pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); const Vector & weights = pair_weight(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); const Vector & to_acc = to_accumulate(first_pair_types->second, ghost_type2); Vector & acc = accumulated(first_pair_types->first, ghost_type1); acc.resize(to_acc.getSize()); acc.clear(); Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::const_iterator< types::Vector > pair_w = weights.begin(2); - typename Vector::template const_iterator< types::Vector > to_acc_it = to_acc.begin(nb_degree_of_freedom); - typename Vector::template iterator< typename types::Vector > acc_it = acc.begin(nb_degree_of_freedom); - for(;first_pair != last_pair; ++first_pair, ++pair_w) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); - for(UInt d = 0; d < nb_degree_of_freedom; ++d) { - acc_it[q1](d) += (*pair_w)(0) * to_acc_it[q2](d); - acc_it[q2](d) += (*pair_w)(1) * to_acc_it[q1](d); + for(UInt d = 0; d < nb_degree_of_freedom; ++d){ + acc(q1, d) += *pair_w * to_acc(q2, d); + acc(q2, d) += *pair_w * to_acc(q1, d); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::createCellList(const ByElementTypeReal & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); const Real safety_factor = 1.2; // for the cell grid spacing Mesh & mesh = this->model->getFEM().getMesh(); mesh.computeBoundingBox(); Real lower_bounds[spatial_dimension]; Real upper_bounds[spatial_dimension]; mesh.getLocalLowerBounds(lower_bounds); mesh.getLocalUpperBounds(upper_bounds); Real spacing[spatial_dimension]; for (UInt i = 0; i < spatial_dimension; ++i) { spacing[i] = radius * safety_factor; } cell_list = new RegularGrid(spatial_dimension, lower_bounds, upper_bounds, spacing); GhostType ghost_type = _not_ghost; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); // first generate the quad points coordinate and count the number of points per cell for(; it != last_type; ++it) { const Vector & quads = quadrature_points_coordinates(*it, ghost_type); Vector::const_iterator first_quad, last_quad; first_quad = quads.begin(spatial_dimension); last_quad = quads.end(spatial_dimension); for(;first_quad != last_quad; ++first_quad) { cell_list->count(*first_quad); } } QuadraturePoint q; q.ghost_type = ghost_type; // second insert the point in the cells cell_list->beginInsertions(); it = mesh.firstType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { Vector & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_quad = this->model->getFEM().getNbQuadraturePoints(*it, ghost_type); const Vector & quads = quadrature_points_coordinates(*it, ghost_type); q.type = *it; Vector::const_iterator quad = quads.begin(spatial_dimension); UInt * elem = elem_filter.storage(); for (UInt e = 0; e < nb_element; ++e) { q.element = *elem; for (UInt nq = 0; nq < nb_quad; ++nq) { q.num_point = nq; q.setPosition(*quad); cell_list->insert(q, *quad); ++quad; } ++elem; } } cell_list->endInsertions(); SynchronizerRegistry & synch_registry = this->model->getSynchronizerRegistry(); std::stringstream sstr; sstr << id << ":grid_synchronizer"; GridSynchronizer * synch = GridSynchronizer::createGridSynchronizer(mesh, *cell_list, sstr.str()); synch_registry.registerSynchronizer(*synch, _gst_mnl_damage); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::updatePairList(const ByElementTypeReal & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); Mesh & mesh = this->model->getFEM().getMesh(); GhostType ghost_type = _not_ghost; // generate the pair of neighbor depending of the cell_list Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { // Preparing datas const Vector & quads = quadrature_points_coordinates(*it, ghost_type); Vector::const_iterator first_quad = quads.begin(spatial_dimension); Vector::const_iterator last_quad = quads.end(spatial_dimension); ByElementTypeUInt & pairs = pair_list(ByElementTypeUInt("pairs", id, memory_id), *it, ghost_type); ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; //Real * neigh_quad_positions = NULL; Vector * neighbors = NULL; UInt my_num_quad = 0; // loop over quad points for(;first_quad != last_quad; ++first_quad, ++my_num_quad) { RegularGrid::Cell cell = cell_list->getCell(*first_quad); RegularGrid::neighbor_cells_iterator first_neigh_cell = cell_list->beginNeighborCells(cell); RegularGrid::neighbor_cells_iterator last_neigh_cell = cell_list->endNeighborCells(cell); // loop over neighbor cells of the one containing the current quadrature // point for (; first_neigh_cell != last_neigh_cell; ++first_neigh_cell) { RegularGrid::iterator first_neigh_quad = cell_list->beginCell(*first_neigh_cell); RegularGrid::iterator last_neigh_quad = cell_list->endCell(*first_neigh_cell); // loop over the quadrature point in the current cell of the cell list for (;first_neigh_quad != last_neigh_quad; ++first_neigh_quad){ QuadraturePoint & quad = *first_neigh_quad; UInt nb_quad_per_elem = this->model->getFEM().getNbQuadraturePoints(quad.type, quad.ghost_type); UInt neigh_num_quad = quad.element * nb_quad_per_elem + quad.num_point; // little optimization to not search in the map at each quad points if(quad.type != current_element_type || quad.ghost_type != current_ghost_type) { // neigh_quad_positions = quadrature_points_coordinates(quad.type, // quad.ghost_type).storage(); current_element_type = quad.type; current_ghost_type = quad.ghost_type; if(!pairs.exists(current_element_type, current_ghost_type)) { neighbors = &(pairs.alloc(0, 2, current_element_type, current_ghost_type)); } else { neighbors = &(pairs(current_element_type, current_ghost_type)); } existing_pairs.insert(std::pair(*it, current_element_type)); } // types::RVector neigh_quad(neigh_quad_positions + neigh_num_quad * spatial_dimension, // spatial_dimension); const types::RVector & neigh_quad = quad.getPosition(); Real distance = first_quad->distance(neigh_quad); if(distance <= radius && my_num_quad <= neigh_num_quad) { // sotring only half lists UInt pair[2]; pair[0] = my_num_quad; pair[1] = neigh_num_quad; neighbors->push_back(pair); } // } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::computeWeights(const ByElementTypeReal & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); std::set< std::pair >::iterator first_pair_types; std::set< std::pair >::iterator last_pair_types = existing_pairs.end(); GhostType ghost_type1, ghost_type2; ghost_type1 = ghost_type2 = _not_ghost; ByElementTypeReal quadrature_points_volumes("quadrature_points_volumes", id, memory_id); this->model->getFEM().getMesh().initByElementTypeVector(quadrature_points_volumes, 1, 0); weight_func->updateInternals(quadrature_points_volumes); // Compute the weights first_pair_types = existing_pairs.begin(); for (; first_pair_types != last_pair_types; ++first_pair_types) { ElementType type1 = first_pair_types->first; ElementType type2 = first_pair_types->second; const Vector & pairs = pair_list(type1, ghost_type1)(type2, ghost_type2); ByElementTypeReal & weights_type_1 = pair_weight(type1, ghost_type1); Vector * tmp_weight = NULL; if(!weights_type_1.exists(type2, ghost_type2)) { tmp_weight = &(weights_type_1.alloc(0, 2, type2, ghost_type2)); } else { tmp_weight = &(weights_type_1(type2, ghost_type2)); } Vector & weights = *tmp_weight; weights.resize(pairs.getSize()); weights.clear(); const Vector & elem_filter = element_filter(type1, ghost_type1); UInt nb_quad1 = this->model->getFEM().getNbQuadraturePoints(type1); UInt nb_quad2 = this->model->getFEM().getNbQuadraturePoints(type2); UInt nb_tot_quad = nb_quad1 * elem_filter.getSize();; Vector & quads_volumes = quadrature_points_volumes(type1, ghost_type1); quads_volumes.resize(nb_tot_quad); quads_volumes.clear(); Vector::const_iterator iquads1; Vector::const_iterator iquads2; iquads1 = quadrature_points_coordinates(type1, ghost_type1).begin(spatial_dimension); iquads2 = quadrature_points_coordinates(type2, ghost_type2).begin(spatial_dimension); Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::iterator weight = weights.begin(2); this->weight_func->selectType(type1, ghost_type1, type2, ghost_type2); // Weight function for(;first_pair != last_pair; ++first_pair, ++weight) { UInt _q1 = (*first_pair)(0); UInt _q2 = (*first_pair)(1); const types::RVector & pos1 = iquads1[_q1]; const types::RVector & pos2 = iquads2[_q2]; QuadraturePoint q1(_q1 / nb_quad1, _q1 % nb_quad1, _q1, pos1, type1, ghost_type1); QuadraturePoint q2(_q2 / nb_quad2, _q2 % nb_quad2, _q2, pos2, type2, ghost_type2); Real r = pos1.distance(pos2); (*weight)(0) = this->weight_func->operator()(r, q1, q2); if(_q1 != _q2) (*weight)(1) = this->weight_func->operator()(r, q2, q1); else (*weight)(1) = 0; quads_volumes(_q1) += (*weight)(0); quads_volumes(_q2) += (*weight)(1); } } //normalize the weights first_pair_types = existing_pairs.begin(); for (; first_pair_types != last_pair_types; ++first_pair_types) { ElementType type1 = first_pair_types->first; ElementType type2 = first_pair_types->second; const Vector & pairs = pair_list(type1, ghost_type1)(type2, ghost_type2); Vector & weights = pair_weight(type1, ghost_type1)(type2, ghost_type2); Vector & quads_volumes = quadrature_points_volumes(type1, ghost_type1); Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::iterator weight = weights.begin(2); for(;first_pair != last_pair; ++first_pair, ++weight) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); (*weight)(0) *= 1. / quads_volumes(q1); (*weight)(1) *= 1. / quads_volumes(q2); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template bool MaterialNonLocal::setParam(const std::string & key, const std::string & value, __attribute__((unused)) const ID & id) { std::stringstream sstr(value); if(key == "radius") { sstr >> radius; } else if(key == "UpdateWeights") { sstr >> update_weigths; } else if(!weight_func->setParam(key, value)) return false; return true; } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::savePairs(const std::string & filename) const { std::ofstream pout; pout.open(filename.c_str()); std::set< std::pair >::const_iterator first_pair_types = existing_pairs.begin(); std::set< std::pair >::const_iterator last_pair_types = existing_pairs.end(); GhostType ghost_type1, ghost_type2; ghost_type1 = ghost_type2 = _not_ghost; for (; first_pair_types != last_pair_types; ++first_pair_types) { const Vector & pairs = pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); const Vector & weights = pair_weight(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); pout << "Types : " << first_pair_types->first << " " << first_pair_types->second << std::endl; Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::const_iterator pair_w = weights.begin(2); for(;first_pair != last_pair; ++first_pair, ++pair_w) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); pout << q1 << " " << q2 << " "<< (*pair_w)(0) << " " << (*pair_w)(1) << std::endl; } } } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::neighbourhoodStatistics(const std::string & filename) const { std::ofstream pout; pout.open(filename.c_str()); std::set< std::pair >::const_iterator first_pair_types = existing_pairs.begin(); std::set< std::pair >::const_iterator last_pair_types = existing_pairs.end(); const Mesh & mesh = this->model->getFEM().getMesh(); ByElementTypeUInt nb_neighbors("nb_neighbours", id, memory_id); initInternalVector(nb_neighbors, 1); resizeInternalVector(nb_neighbors); GhostType ghost_type1, ghost_type2; ghost_type1 = ghost_type2 = _not_ghost; for (; first_pair_types != last_pair_types; ++first_pair_types) { const Vector & pairs = pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); pout << "Types : " << first_pair_types->first << " " << first_pair_types->second << std::endl; Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector & nb_neigh_1 = nb_neighbors(first_pair_types->first, ghost_type1); Vector & nb_neigh_2 = nb_neighbors(first_pair_types->second, ghost_type2); for(;first_pair != last_pair; ++first_pair) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); ++(nb_neigh_1(q1)); if(q1 != q2) ++(nb_neigh_2(q2)); } } Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type1); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type1); UInt nb_quads = 0; Real mean_nb_neig = 0; UInt max_nb_neig = 0; UInt min_nb_neig = std::numeric_limits::max(); for(; it != last_type; ++it) { Vector & nb_neighor = nb_neighbors(*it, ghost_type1); Vector::iterator nb_neigh = nb_neighor.begin(); Vector::iterator end_neigh = nb_neighor.end(); for (; nb_neigh != end_neigh; ++nb_neigh, ++nb_quads) { UInt nb = *nb_neigh; mean_nb_neig += nb; max_nb_neig = std::max(max_nb_neig, nb); min_nb_neig = std::min(min_nb_neig, nb); } } mean_nb_neig /= Real(nb_quads); pout << "Nb quadrature points: " << nb_quads << std::endl; pout << "Average nb neighbors: " << mean_nb_neig << std::endl; pout << "Max nb neighbors: " << max_nb_neig << std::endl; pout << "Min nb neighbors: " << min_nb_neig << std::endl; pout.close(); } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Material<_non_local> [" << std::endl; stream << space << " + Radius : " << radius << std::endl; stream << space << " + UpdateWeights : " << update_weigths << std::endl; stream << space << " + Weight Function : " << *weight_func << std::endl; stream << space << "]" << std::endl; } diff --git a/test/test_model/test_solid_mechanics_model/patch_tests/implicit/patch_test_implicit.cc b/test/test_model/test_solid_mechanics_model/patch_tests/implicit/patch_test_implicit.cc index 035a0d165..1455a8bd9 100644 --- a/test/test_model/test_solid_mechanics_model/patch_tests/implicit/patch_test_implicit.cc +++ b/test/test_model/test_solid_mechanics_model/patch_tests/implicit/patch_test_implicit.cc @@ -1,308 +1,308 @@ /** * @file test_check_stress.cc * @author David Kammer * @date Wed Feb 16 13:56:42 2011 * * @brief patch test for elastic material in solid mechanics model * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "element_class.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.hh" #endif //AKANTU_USE_IOHELPER using namespace akantu; Real alpha [3][4] = { { 0.01, 0.02, 0.03, 0.04 }, { 0.05, 0.06, 0.07, 0.08 }, { 0.09, 0.10, 0.11, 0.12 } }; /* -------------------------------------------------------------------------- */ template static types::Matrix prescribed_strain() { UInt spatial_dimension = ElementClass::getSpatialDimension(); types::Matrix strain(spatial_dimension, spatial_dimension); for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = 0; j < spatial_dimension; ++j) { strain(i, j) = alpha[i][j + 1]; } } return strain; } template static types::Matrix prescribed_stress() { UInt spatial_dimension = ElementClass::getSpatialDimension(); types::Matrix stress(spatial_dimension, spatial_dimension); //plane strain in 2d types::Matrix strain(spatial_dimension, spatial_dimension); types::Matrix pstrain; pstrain = prescribed_strain(); Real nu = 0.3; Real E = 2.1e11; Real trace = 0; /// symetric part of the strain tensor for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) strain(i,j) = 0.5 * (pstrain(i, j) + pstrain(j, i)); for (UInt i = 0; i < spatial_dimension; ++i) trace += strain(i,i); Real Ep = E / (1 + nu); for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) { stress(i, j) = Ep * strain(i,j); if(i == j) stress(i, j) += Ep * (nu / (1 - 2*nu)) * trace; } return stress; } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER static void paraviewInit(iohelper::Dumper & dumper, const SolidMechanicsModel & model); static void paraviewDump(iohelper::Dumper & dumper); #endif /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { initialize(argc, argv); UInt dim = ElementClass::getSpatialDimension(); const ElementType element_type = TYPE; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; std::stringstream filename; filename << TYPE << ".msh"; mesh_io.read(filename.str(), my_mesh); UInt nb_nodes = my_mesh.getNbNodes(); /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors my_model.getForce().clear(); my_model.getVelocity().clear(); my_model.getAcceleration().clear(); my_model.getDisplacement().clear(); my_model.initModel(); my_model.readMaterials("material_check_stress.dat"); my_model.initMaterials(); my_model.initImplicit(); const Vector & coordinates = my_mesh.getNodes(); Vector & displacement = my_model.getDisplacement(); Vector & boundary = my_model.getBoundary(); MeshUtils::buildFacets(my_mesh); MeshUtils::buildSurfaceID(my_mesh); CSR surface_nodes; MeshUtils::buildNodesPerSurface(my_mesh, surface_nodes); for (UInt s = 0; s < surface_nodes.getNbRows(); ++s) { CSR::iterator snode = surface_nodes.begin(s); for(; snode != surface_nodes.end(s); ++snode) { UInt n = *snode; std::cout << "Node " << n << std::endl; for (UInt i = 0; i < dim; ++i) { displacement(n, i) = alpha[i][0]; for (UInt j = 0; j < dim; ++j) { displacement(n, i) += alpha[i][j + 1] * coordinates(n, j); } boundary(n, i) = true; } } } /* ------------------------------------------------------------------------ */ /* Main loop */ /* ------------------------------------------------------------------------ */ akantu::UInt count = 0; my_model.updateResidual(); #ifdef AKANTU_USE_IOHELPER iohelper::DumperParaview dumper; paraviewInit(dumper, my_model); #endif while(!my_model.testConvergenceResidual(2e-4) && (count < 100)) { std::cout << "Iter : " << ++count << std::endl; my_model.assembleStiffnessMatrix(); my_model.solveStatic(); my_model.updateResidual(); } #ifdef AKANTU_USE_IOHELPER paraviewDump(dumper); #endif if(count > 1) { std::cerr << "The code did not converge in 1 step !" << std::endl; return EXIT_FAILURE; } /* ------------------------------------------------------------------------ */ /* Checks */ /* ------------------------------------------------------------------------ */ UInt nb_quadrature_points = my_model.getFEM().getNbQuadraturePoints(element_type); Vector & stress_vect = const_cast &>(my_model.getMaterial(0).getStress(element_type)); Vector & strain_vect = const_cast &>(my_model.getMaterial(0).getStrain(element_type)); Vector::iterator stress_it = stress_vect.begin(dim, dim); Vector::iterator strain_it = strain_vect.begin(dim, dim); types::Matrix presc_stress; presc_stress = prescribed_stress(); types::Matrix presc_strain; presc_strain = prescribed_strain(); UInt nb_element = my_mesh.getNbElement(TYPE); for (UInt el = 0; el < nb_element; ++el) { for (UInt q = 0; q < nb_quadrature_points; ++q) { types::Matrix & stress = *stress_it; types::Matrix & strain = *strain_it; for (UInt i = 0; i < dim; ++i) { for (UInt j = 0; j < dim; ++j) { if(!(std::abs(strain(i, j) - presc_strain(i, j)) < 1e-15)) { std::cerr << "strain[" << i << "," << j << "] = " << strain(i, j) << " but should be = " << presc_strain(i, j) << " (-" << std::abs(strain(i, j) - presc_strain(i, j)) << ") [el : " << el<< " - q : " << q << "]" << std::endl; std::cerr << "computed : " << strain << "reference : " << presc_strain << std::endl; return EXIT_FAILURE; } if(!(std::abs(stress(i, j) - presc_stress(i, j)) < 1e-3)) { std::cerr << "stress[" << i << "," << j << "] = " << stress(i, j) << " but should be = " << presc_stress(i, j) << " (-" << std::abs(stress(i, j) - presc_stress(i, j)) << ") [el : " << el<< " - q : " << q << "]" << std::endl; std::cerr << "computed : " << stress << "reference : " << presc_stress << std::endl; return EXIT_FAILURE; } } } ++stress_it; ++strain_it; } } for (UInt n = 0; n < nb_nodes; ++n) { for (UInt i = 0; i < dim; ++i) { Real disp = alpha[i][0]; for (UInt j = 0; j < dim; ++j) { disp += alpha[i][j + 1] * coordinates(n, j); } if(!(std::abs(displacement(n,i) - disp) < 1e-15)) { std::cerr << "displacement(" << n << ", " << i <<")=" << displacement(n,i) << " should be equal to " << disp << std::endl; return EXIT_FAILURE; } } } // std::cout << "Strain : " << strain; // std::cout << "Stress : " << stress; // finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ /* iohelper::Dumper vars */ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER template static iohelper::ElemType paraviewType(); -template <> static iohelper::ElemType paraviewType<_segment_2>() { return iohelper::LINE1; } -template <> static iohelper::ElemType paraviewType<_segment_3>() { return iohelper::LINE2; } -template <> static iohelper::ElemType paraviewType<_triangle_3>() { return iohelper::TRIANGLE1; } -template <> static iohelper::ElemType paraviewType<_triangle_6>() { return iohelper::TRIANGLE2; } -template <> static iohelper::ElemType paraviewType<_quadrangle_4>() { return iohelper::QUAD1; } -template <> static iohelper::ElemType paraviewType<_quadrangle_8>() { return iohelper::QUAD2; } -template <> static iohelper::ElemType paraviewType<_tetrahedron_4>() { return iohelper::TETRA1; } -template <> static iohelper::ElemType paraviewType<_tetrahedron_10>() { return iohelper::TETRA2; } -template <> static iohelper::ElemType paraviewType<_hexahedron_8>() { return iohelper::HEX1; } +template <> iohelper::ElemType paraviewType<_segment_2>() { return iohelper::LINE1; } +template <> iohelper::ElemType paraviewType<_segment_3>() { return iohelper::LINE2; } +template <> iohelper::ElemType paraviewType<_triangle_3>() { return iohelper::TRIANGLE1; } +template <> iohelper::ElemType paraviewType<_triangle_6>() { return iohelper::TRIANGLE2; } +template <> iohelper::ElemType paraviewType<_quadrangle_4>() { return iohelper::QUAD1; } +template <> iohelper::ElemType paraviewType<_quadrangle_8>() { return iohelper::QUAD2; } +template <> iohelper::ElemType paraviewType<_tetrahedron_4>() { return iohelper::TETRA1; } +template <> iohelper::ElemType paraviewType<_tetrahedron_10>() { return iohelper::TETRA2; } +template <> iohelper::ElemType paraviewType<_hexahedron_8>() { return iohelper::HEX1; } void paraviewInit(iohelper::Dumper & dumper, const SolidMechanicsModel & model) { UInt spatial_dimension = ElementClass::getSpatialDimension(); UInt nb_nodes = model.getFEM().getMesh().getNbNodes(); UInt nb_element = model.getFEM().getMesh().getNbElement(TYPE); std::stringstream filename; filename << "out_" << TYPE; dumper.SetMode(iohelper::TEXT); dumper.SetPoints(model.getFEM().getMesh().getNodes().values, spatial_dimension, nb_nodes, filename.str().c_str()); dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(TYPE).values, paraviewType(), nb_element, iohelper::C_MODE); dumper.AddNodeDataField(model.getDisplacement().values, spatial_dimension, "displacements"); dumper.AddNodeDataField(model.getVelocity().values, spatial_dimension, "velocity"); dumper.AddNodeDataField(model.getResidual().values, spatial_dimension, "force"); dumper.AddNodeDataField(model.getMass().values, spatial_dimension, "mass"); dumper.AddNodeDataField(model.getForce().values, spatial_dimension, "applied_force"); dumper.AddElemDataField(model.getMaterial(0).getStrain(TYPE).values, spatial_dimension*spatial_dimension, "strain"); dumper.AddElemDataField(model.getMaterial(0).getStrain(TYPE).values, spatial_dimension*spatial_dimension, "stress"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("applied_force", 1); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); } /* -------------------------------------------------------------------------- */ void paraviewDump(iohelper::Dumper & dumper) { dumper.Dump(); } #endif diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/io_helper_tools.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/io_helper_tools.cc index 548f9ee05..230775a33 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/io_helper_tools.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/io_helper_tools.cc @@ -1,294 +1,294 @@ /** * @file io_helper_tools.cc * @author Nicolas Richart * @date Fri Sep 30 11:13:01 2011 * * @brief * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "io_helper_tools.hh" #include "aka_common.hh" #include "mesh.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "static_communicator.hh" // #include "dof_synchronizer.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; /* ------------------------------------------------------------------------ */ template static iohelper::ElemType getIOHelperType() { AKANTU_DEBUG_TO_IMPLEMENT(); return iohelper::MAX_ELEM_TYPE; } -template <> static iohelper::ElemType getIOHelperType<_segment_2>() { return iohelper::LINE1; } -template <> static iohelper::ElemType getIOHelperType<_segment_3>() { return iohelper::LINE2; } -template <> static iohelper::ElemType getIOHelperType<_triangle_3>() { return iohelper::TRIANGLE1; } -template <> static iohelper::ElemType getIOHelperType<_triangle_6>() { return iohelper::TRIANGLE2; } -template <> static iohelper::ElemType getIOHelperType<_quadrangle_4>() { return iohelper::QUAD1; } -template <> static iohelper::ElemType getIOHelperType<_quadrangle_8>() { return iohelper::QUAD2; } -template <> static iohelper::ElemType getIOHelperType<_tetrahedron_4>() { return iohelper::TETRA1; } -template <> static iohelper::ElemType getIOHelperType<_tetrahedron_10>() { return iohelper::TETRA2; } -template <> static iohelper::ElemType getIOHelperType<_hexahedron_8>() { return iohelper::HEX1; } +template <> iohelper::ElemType getIOHelperType<_segment_2>() { return iohelper::LINE1; } +template <> iohelper::ElemType getIOHelperType<_segment_3>() { return iohelper::LINE2; } +template <> iohelper::ElemType getIOHelperType<_triangle_3>() { return iohelper::TRIANGLE1; } +template <> iohelper::ElemType getIOHelperType<_triangle_6>() { return iohelper::TRIANGLE2; } +template <> iohelper::ElemType getIOHelperType<_quadrangle_4>() { return iohelper::QUAD1; } +template <> iohelper::ElemType getIOHelperType<_quadrangle_8>() { return iohelper::QUAD2; } +template <> iohelper::ElemType getIOHelperType<_tetrahedron_4>() { return iohelper::TETRA1; } +template <> iohelper::ElemType getIOHelperType<_tetrahedron_10>() { return iohelper::TETRA2; } +template <> iohelper::ElemType getIOHelperType<_hexahedron_8>() { return iohelper::HEX1; } static iohelper::ElemType getIOHelperType(ElementType type) { iohelper::ElemType ioh_type; #define GET_IOHELPER_TYPE(type) \ ioh_type = getIOHelperType(); AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(GET_IOHELPER_TYPE); #undef GET_IOHELPER_TYPE return ioh_type; } /* ------------------------------------------------------------------------ */ void paraviewInit(iohelper::Dumper & dumper, const SolidMechanicsModel & model, const ElementType & type, const std::string & filename) { const Mesh & mesh = model.getFEM().getMesh(); UInt spatial_dimension = mesh.getSpatialDimension(type); UInt nb_nodes = mesh.getNbNodes(); UInt nb_element = mesh.getNbElement(type); std::stringstream filename_sstr; filename_sstr << filename << "_" << type; UInt whoami = StaticCommunicator::getStaticCommunicator()->whoAmI(); UInt nproc = StaticCommunicator::getStaticCommunicator()->getNbProc(); dumper.SetMode(iohelper::TEXT); dumper.SetParallelContext(whoami, nproc); dumper.SetPoints(mesh.getNodes().values, spatial_dimension, nb_nodes, filename_sstr.str().c_str()); dumper.SetConnectivity((int *)mesh.getConnectivity(type).values, getIOHelperType(type), nb_element, iohelper::C_MODE); dumper.AddNodeDataField(model.getDisplacement().values, spatial_dimension, "displacements"); dumper.AddNodeDataField(model.getVelocity().values, spatial_dimension, "velocity"); dumper.AddNodeDataField(model.getAcceleration().values, spatial_dimension, "acceleration"); dumper.AddNodeDataField(model.getMass().values, spatial_dimension, "mass"); dumper.AddNodeDataField(model.getResidual().values, spatial_dimension, "force"); dumper.AddNodeDataField(model.getForce().values, spatial_dimension, "applied_force"); dumper.AddElemDataField(model.getMaterial(0).getStrain(type).values, spatial_dimension*spatial_dimension, "strain"); dumper.AddElemDataField(model.getMaterial(0).getStress(type).values, spatial_dimension*spatial_dimension, "stress"); if(dynamic_cast(&model.getMaterial(0))) { const MaterialDamage & mat = dynamic_cast(model.getMaterial(0)); Real * dam = mat.getDamage(type).storage(); dumper.AddElemDataField(dam, 1, "damage"); } dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("applied_force", 1); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); } /* ------------------------------------------------------------------------ */ void paraviewDump(iohelper::Dumper & dumper) { dumper.Dump(); } /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ // Vector checkpoint_displacements(0, spatial_dimension); // Vector checkpoint_velocity(0, spatial_dimension); // Vector checkpoint_acceleration(0, spatial_dimension); // Vector checkpoint_force(0, spatial_dimension); // /* ------------------------------------------------------------------------ */ // void checkpointInit(iohelper::Dumper & dumper, // const SolidMechanicsModel & model, // const ElementType & type, // const std::string & filename) { // UInt whoami = StaticCommunicator::getStaticCommunicator()->whoAmI(); // UInt nproc = StaticCommunicator::getStaticCommunicator()->getNbProc(); // Vector & displacements = model.getDisplacement(); // Vector & velocity = model.getVelocity(); // Vector & acceleration = model.getAcceleration(); // Vector & force = model.getForce(); // DOFSynchronizer & dof_synchronizer = const_cast(model.getDOFSynchronizer()); // dof_synchronizer.initScatterGatherCommunicationScheme(); // if(whoami == 0){ // const Mesh & mesh = model.getFEM().getMesh(); // UInt nb_nodes = mesh.getNbGlobalNodes(); // checkpoint_displacements.resize(nb_nodes); // checkpoint_velocity .resize(nb_nodes); // checkpoint_acceleration .resize(nb_nodes); // checkpoint_force .resize(nb_nodes); // dof_synchronizer.gather(displacements, 0, &checkpoint_displacements); // dof_synchronizer.gather(velocity , 0, &checkpoint_velocity ); // dof_synchronizer.gather(acceleration , 0, &checkpoint_acceleration ); // dof_synchronizer.gather(force , 0, &checkpoint_force ); // UInt nb_element = mesh.getNbElement(type); // UInt spatial_dimension = mesh.getSpatialDimension(type); // std::stringstream filename_sstr; filename_sstr << filename << "_" << type; // dumper.SetMode(iohelper::COMPRESSED); // dumper.SetParallelContext(whoami, nproc); // dumper.SetPoints(mesh.getNodes().values, // spatial_dimension, nb_nodes, filename_sstr.str().c_str()); // dumper.SetConnectivity((int *)mesh.getConnectivity(type).values, // getIOHelperType(type), nb_element, iohelper::C_MODE); // dumper.AddNodeDataField(checkpoint_displacements.storage(), spatial_dimension, "displacements"); // dumper.AddNodeDataField(checkpoint_velocity .storage(), spatial_dimension, "velocity"); // dumper.AddNodeDataField(checkpoint_acceleration .storage(), spatial_dimension, "acceleration"); // dumper.AddNodeDataField(checkpoint_force .storage(), spatial_dimension, "applied_force"); // dumper.SetPrefix("restart/"); // dumper.Init(); // dumper.Dump(); // } else { // dof_synchronizer.gather(displacements, 0); // dof_synchronizer.gather(velocity , 0); // dof_synchronizer.gather(acceleration , 0); // dof_synchronizer.gather(force , 0); // } // } // /* ------------------------------------------------------------------------ */ // void checkpoint(iohelper::Dumper & dumper, // const SolidMechanicsModel & model) { // UInt whoami = StaticCommunicator::getStaticCommunicator()->whoAmI(); // DOFSynchronizer & dof_synchronizer = const_cast(model.getDOFSynchronizer()); // Vector & displacements = model.getDisplacement(); // Vector & velocity = model.getVelocity(); // Vector & acceleration = model.getAcceleration(); // Vector & force = model.getForce(); // if(whoami == 0){ // dof_synchronizer.gather(displacements, 0, &checkpoint_displacements); // dof_synchronizer.gather(velocity , 0, &checkpoint_velocity ); // dof_synchronizer.gather(acceleration , 0, &checkpoint_acceleration ); // dof_synchronizer.gather(force , 0, &checkpoint_force ); // dumper.Dump(); // } else { // dof_synchronizer.gather(displacements, 0); // dof_synchronizer.gather(velocity , 0); // dof_synchronizer.gather(acceleration , 0); // dof_synchronizer.gather(force , 0); // } // } // /* ------------------------------------------------------------------------ */ // /* ------------------------------------------------------------------------ */ // void restart(const SolidMechanicsModel & model, // const ElementType & type, // const std::string & filename) { // UInt whoami = StaticCommunicator::getStaticCommunicator()->whoAmI(); // UInt nproc = StaticCommunicator::getStaticCommunicator()->getNbProc(); // DOFSynchronizer & dof_synchronizer = const_cast(model.getDOFSynchronizer()); // dof_synchronizer.initScatterGatherCommunicationScheme(); // Vector & displacements = model.getDisplacement(); // Vector & velocity = model.getVelocity(); // Vector & acceleration = model.getAcceleration(); // Vector & force = model.getForce(); // if(whoami == 0){ // const Mesh & mesh = model.getFEM().getMesh(); // UInt nb_nodes = mesh.getNbGlobalNodes(); // UInt spatial_dimension = mesh.getSpatialDimension(type); // std::stringstream filename_sstr; filename_sstr << filename << "_" << type; // iohelper::ReaderRestart reader; // reader.SetMode(iohelper::COMPRESSED); // reader.SetParallelContext(whoami, nproc); // reader.SetPoints(filename_sstr.str().c_str()); // reader.SetConnectivity(getIOHelperType(type)); // reader.AddNodeDataField("displacements"); // reader.AddNodeDataField("velocity"); // reader.AddNodeDataField("acceleration"); // reader.AddNodeDataField("applied_force"); // reader.SetPrefix("restart/"); // reader.Init(); // reader.Read(); // Vector restart_tmp_vect(nb_nodes, spatial_dimension); // displacements.clear(); // velocity.clear(); // acceleration.clear(); // force.clear(); // Real * tmp = reader.GetNodeDataField("displacements"); // std::copy(tmp, tmp + nb_nodes * spatial_dimension, restart_tmp_vect.storage()); // dof_synchronizer.scatter(displacements, 0, &restart_tmp_vect); // tmp = reader.GetNodeDataField("velocity"); // std::copy(tmp, tmp + nb_nodes * spatial_dimension, restart_tmp_vect.storage()); // dof_synchronizer.scatter(velocity , 0, &restart_tmp_vect); // tmp = reader.GetNodeDataField("acceleration"); // std::copy(tmp, tmp + nb_nodes * spatial_dimension, restart_tmp_vect.storage()); // dof_synchronizer.scatter(acceleration , 0, &restart_tmp_vect); // tmp = reader.GetNodeDataField("applied_force"); // std::copy(tmp, tmp + nb_nodes * spatial_dimension, restart_tmp_vect.storage()); // dof_synchronizer.scatter(force , 0, &restart_tmp_vect); // } else { // displacements.clear(); // velocity.clear(); // acceleration.clear(); // force.clear(); // dof_synchronizer.scatter(displacements, 0); // dof_synchronizer.scatter(velocity , 0); // dof_synchronizer.scatter(acceleration , 0); // dof_synchronizer.scatter(force , 0); // } // } /* ------------------------------------------------------------------------ */ diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/test_contact_search_explicit/io_helper_tools.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/test_contact_search_explicit/io_helper_tools.cc index 548f9ee05..230775a33 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/test_contact_search_explicit/io_helper_tools.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/test_contact_search_explicit/io_helper_tools.cc @@ -1,294 +1,294 @@ /** * @file io_helper_tools.cc * @author Nicolas Richart * @date Fri Sep 30 11:13:01 2011 * * @brief * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "io_helper_tools.hh" #include "aka_common.hh" #include "mesh.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "static_communicator.hh" // #include "dof_synchronizer.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; /* ------------------------------------------------------------------------ */ template static iohelper::ElemType getIOHelperType() { AKANTU_DEBUG_TO_IMPLEMENT(); return iohelper::MAX_ELEM_TYPE; } -template <> static iohelper::ElemType getIOHelperType<_segment_2>() { return iohelper::LINE1; } -template <> static iohelper::ElemType getIOHelperType<_segment_3>() { return iohelper::LINE2; } -template <> static iohelper::ElemType getIOHelperType<_triangle_3>() { return iohelper::TRIANGLE1; } -template <> static iohelper::ElemType getIOHelperType<_triangle_6>() { return iohelper::TRIANGLE2; } -template <> static iohelper::ElemType getIOHelperType<_quadrangle_4>() { return iohelper::QUAD1; } -template <> static iohelper::ElemType getIOHelperType<_quadrangle_8>() { return iohelper::QUAD2; } -template <> static iohelper::ElemType getIOHelperType<_tetrahedron_4>() { return iohelper::TETRA1; } -template <> static iohelper::ElemType getIOHelperType<_tetrahedron_10>() { return iohelper::TETRA2; } -template <> static iohelper::ElemType getIOHelperType<_hexahedron_8>() { return iohelper::HEX1; } +template <> iohelper::ElemType getIOHelperType<_segment_2>() { return iohelper::LINE1; } +template <> iohelper::ElemType getIOHelperType<_segment_3>() { return iohelper::LINE2; } +template <> iohelper::ElemType getIOHelperType<_triangle_3>() { return iohelper::TRIANGLE1; } +template <> iohelper::ElemType getIOHelperType<_triangle_6>() { return iohelper::TRIANGLE2; } +template <> iohelper::ElemType getIOHelperType<_quadrangle_4>() { return iohelper::QUAD1; } +template <> iohelper::ElemType getIOHelperType<_quadrangle_8>() { return iohelper::QUAD2; } +template <> iohelper::ElemType getIOHelperType<_tetrahedron_4>() { return iohelper::TETRA1; } +template <> iohelper::ElemType getIOHelperType<_tetrahedron_10>() { return iohelper::TETRA2; } +template <> iohelper::ElemType getIOHelperType<_hexahedron_8>() { return iohelper::HEX1; } static iohelper::ElemType getIOHelperType(ElementType type) { iohelper::ElemType ioh_type; #define GET_IOHELPER_TYPE(type) \ ioh_type = getIOHelperType(); AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(GET_IOHELPER_TYPE); #undef GET_IOHELPER_TYPE return ioh_type; } /* ------------------------------------------------------------------------ */ void paraviewInit(iohelper::Dumper & dumper, const SolidMechanicsModel & model, const ElementType & type, const std::string & filename) { const Mesh & mesh = model.getFEM().getMesh(); UInt spatial_dimension = mesh.getSpatialDimension(type); UInt nb_nodes = mesh.getNbNodes(); UInt nb_element = mesh.getNbElement(type); std::stringstream filename_sstr; filename_sstr << filename << "_" << type; UInt whoami = StaticCommunicator::getStaticCommunicator()->whoAmI(); UInt nproc = StaticCommunicator::getStaticCommunicator()->getNbProc(); dumper.SetMode(iohelper::TEXT); dumper.SetParallelContext(whoami, nproc); dumper.SetPoints(mesh.getNodes().values, spatial_dimension, nb_nodes, filename_sstr.str().c_str()); dumper.SetConnectivity((int *)mesh.getConnectivity(type).values, getIOHelperType(type), nb_element, iohelper::C_MODE); dumper.AddNodeDataField(model.getDisplacement().values, spatial_dimension, "displacements"); dumper.AddNodeDataField(model.getVelocity().values, spatial_dimension, "velocity"); dumper.AddNodeDataField(model.getAcceleration().values, spatial_dimension, "acceleration"); dumper.AddNodeDataField(model.getMass().values, spatial_dimension, "mass"); dumper.AddNodeDataField(model.getResidual().values, spatial_dimension, "force"); dumper.AddNodeDataField(model.getForce().values, spatial_dimension, "applied_force"); dumper.AddElemDataField(model.getMaterial(0).getStrain(type).values, spatial_dimension*spatial_dimension, "strain"); dumper.AddElemDataField(model.getMaterial(0).getStress(type).values, spatial_dimension*spatial_dimension, "stress"); if(dynamic_cast(&model.getMaterial(0))) { const MaterialDamage & mat = dynamic_cast(model.getMaterial(0)); Real * dam = mat.getDamage(type).storage(); dumper.AddElemDataField(dam, 1, "damage"); } dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("applied_force", 1); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); } /* ------------------------------------------------------------------------ */ void paraviewDump(iohelper::Dumper & dumper) { dumper.Dump(); } /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ // Vector checkpoint_displacements(0, spatial_dimension); // Vector checkpoint_velocity(0, spatial_dimension); // Vector checkpoint_acceleration(0, spatial_dimension); // Vector checkpoint_force(0, spatial_dimension); // /* ------------------------------------------------------------------------ */ // void checkpointInit(iohelper::Dumper & dumper, // const SolidMechanicsModel & model, // const ElementType & type, // const std::string & filename) { // UInt whoami = StaticCommunicator::getStaticCommunicator()->whoAmI(); // UInt nproc = StaticCommunicator::getStaticCommunicator()->getNbProc(); // Vector & displacements = model.getDisplacement(); // Vector & velocity = model.getVelocity(); // Vector & acceleration = model.getAcceleration(); // Vector & force = model.getForce(); // DOFSynchronizer & dof_synchronizer = const_cast(model.getDOFSynchronizer()); // dof_synchronizer.initScatterGatherCommunicationScheme(); // if(whoami == 0){ // const Mesh & mesh = model.getFEM().getMesh(); // UInt nb_nodes = mesh.getNbGlobalNodes(); // checkpoint_displacements.resize(nb_nodes); // checkpoint_velocity .resize(nb_nodes); // checkpoint_acceleration .resize(nb_nodes); // checkpoint_force .resize(nb_nodes); // dof_synchronizer.gather(displacements, 0, &checkpoint_displacements); // dof_synchronizer.gather(velocity , 0, &checkpoint_velocity ); // dof_synchronizer.gather(acceleration , 0, &checkpoint_acceleration ); // dof_synchronizer.gather(force , 0, &checkpoint_force ); // UInt nb_element = mesh.getNbElement(type); // UInt spatial_dimension = mesh.getSpatialDimension(type); // std::stringstream filename_sstr; filename_sstr << filename << "_" << type; // dumper.SetMode(iohelper::COMPRESSED); // dumper.SetParallelContext(whoami, nproc); // dumper.SetPoints(mesh.getNodes().values, // spatial_dimension, nb_nodes, filename_sstr.str().c_str()); // dumper.SetConnectivity((int *)mesh.getConnectivity(type).values, // getIOHelperType(type), nb_element, iohelper::C_MODE); // dumper.AddNodeDataField(checkpoint_displacements.storage(), spatial_dimension, "displacements"); // dumper.AddNodeDataField(checkpoint_velocity .storage(), spatial_dimension, "velocity"); // dumper.AddNodeDataField(checkpoint_acceleration .storage(), spatial_dimension, "acceleration"); // dumper.AddNodeDataField(checkpoint_force .storage(), spatial_dimension, "applied_force"); // dumper.SetPrefix("restart/"); // dumper.Init(); // dumper.Dump(); // } else { // dof_synchronizer.gather(displacements, 0); // dof_synchronizer.gather(velocity , 0); // dof_synchronizer.gather(acceleration , 0); // dof_synchronizer.gather(force , 0); // } // } // /* ------------------------------------------------------------------------ */ // void checkpoint(iohelper::Dumper & dumper, // const SolidMechanicsModel & model) { // UInt whoami = StaticCommunicator::getStaticCommunicator()->whoAmI(); // DOFSynchronizer & dof_synchronizer = const_cast(model.getDOFSynchronizer()); // Vector & displacements = model.getDisplacement(); // Vector & velocity = model.getVelocity(); // Vector & acceleration = model.getAcceleration(); // Vector & force = model.getForce(); // if(whoami == 0){ // dof_synchronizer.gather(displacements, 0, &checkpoint_displacements); // dof_synchronizer.gather(velocity , 0, &checkpoint_velocity ); // dof_synchronizer.gather(acceleration , 0, &checkpoint_acceleration ); // dof_synchronizer.gather(force , 0, &checkpoint_force ); // dumper.Dump(); // } else { // dof_synchronizer.gather(displacements, 0); // dof_synchronizer.gather(velocity , 0); // dof_synchronizer.gather(acceleration , 0); // dof_synchronizer.gather(force , 0); // } // } // /* ------------------------------------------------------------------------ */ // /* ------------------------------------------------------------------------ */ // void restart(const SolidMechanicsModel & model, // const ElementType & type, // const std::string & filename) { // UInt whoami = StaticCommunicator::getStaticCommunicator()->whoAmI(); // UInt nproc = StaticCommunicator::getStaticCommunicator()->getNbProc(); // DOFSynchronizer & dof_synchronizer = const_cast(model.getDOFSynchronizer()); // dof_synchronizer.initScatterGatherCommunicationScheme(); // Vector & displacements = model.getDisplacement(); // Vector & velocity = model.getVelocity(); // Vector & acceleration = model.getAcceleration(); // Vector & force = model.getForce(); // if(whoami == 0){ // const Mesh & mesh = model.getFEM().getMesh(); // UInt nb_nodes = mesh.getNbGlobalNodes(); // UInt spatial_dimension = mesh.getSpatialDimension(type); // std::stringstream filename_sstr; filename_sstr << filename << "_" << type; // iohelper::ReaderRestart reader; // reader.SetMode(iohelper::COMPRESSED); // reader.SetParallelContext(whoami, nproc); // reader.SetPoints(filename_sstr.str().c_str()); // reader.SetConnectivity(getIOHelperType(type)); // reader.AddNodeDataField("displacements"); // reader.AddNodeDataField("velocity"); // reader.AddNodeDataField("acceleration"); // reader.AddNodeDataField("applied_force"); // reader.SetPrefix("restart/"); // reader.Init(); // reader.Read(); // Vector restart_tmp_vect(nb_nodes, spatial_dimension); // displacements.clear(); // velocity.clear(); // acceleration.clear(); // force.clear(); // Real * tmp = reader.GetNodeDataField("displacements"); // std::copy(tmp, tmp + nb_nodes * spatial_dimension, restart_tmp_vect.storage()); // dof_synchronizer.scatter(displacements, 0, &restart_tmp_vect); // tmp = reader.GetNodeDataField("velocity"); // std::copy(tmp, tmp + nb_nodes * spatial_dimension, restart_tmp_vect.storage()); // dof_synchronizer.scatter(velocity , 0, &restart_tmp_vect); // tmp = reader.GetNodeDataField("acceleration"); // std::copy(tmp, tmp + nb_nodes * spatial_dimension, restart_tmp_vect.storage()); // dof_synchronizer.scatter(acceleration , 0, &restart_tmp_vect); // tmp = reader.GetNodeDataField("applied_force"); // std::copy(tmp, tmp + nb_nodes * spatial_dimension, restart_tmp_vect.storage()); // dof_synchronizer.scatter(force , 0, &restart_tmp_vect); // } else { // displacements.clear(); // velocity.clear(); // acceleration.clear(); // force.clear(); // dof_synchronizer.scatter(displacements, 0); // dof_synchronizer.scatter(velocity , 0); // dof_synchronizer.scatter(acceleration , 0); // dof_synchronizer.scatter(force , 0); // } // } /* ------------------------------------------------------------------------ */ diff --git a/test/test_model/test_solid_mechanics_model/test_materials/io_helper_tools.cc b/test/test_model/test_solid_mechanics_model/test_materials/io_helper_tools.cc index 5062962f2..5c4f78c57 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/io_helper_tools.cc +++ b/test/test_model/test_solid_mechanics_model/test_materials/io_helper_tools.cc @@ -1,294 +1,294 @@ /** * @file io_helper_tools.cc * @author Nicolas Richart * @date Fri Sep 30 11:13:01 2011 * * @brief * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "io_helper_tools.hh" #include "aka_common.hh" #include "mesh.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "static_communicator.hh" // #include "dof_synchronizer.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; /* ------------------------------------------------------------------------ */ template static iohelper::ElemType getIOHelperType() { AKANTU_DEBUG_TO_IMPLEMENT(); return iohelper::MAX_ELEM_TYPE; } -template <> static iohelper::ElemType getIOHelperType<_segment_2>() { return iohelper::LINE1; } -template <> static iohelper::ElemType getIOHelperType<_segment_3>() { return iohelper::LINE2; } -template <> static iohelper::ElemType getIOHelperType<_triangle_3>() { return iohelper::TRIANGLE1; } -template <> static iohelper::ElemType getIOHelperType<_triangle_6>() { return iohelper::TRIANGLE2; } -template <> static iohelper::ElemType getIOHelperType<_quadrangle_4>() { return iohelper::QUAD1; } -template <> static iohelper::ElemType getIOHelperType<_quadrangle_8>() { return iohelper::QUAD2; } -template <> static iohelper::ElemType getIOHelperType<_tetrahedron_4>() { return iohelper::TETRA1; } -template <> static iohelper::ElemType getIOHelperType<_tetrahedron_10>() { return iohelper::TETRA2; } -template <> static iohelper::ElemType getIOHelperType<_hexahedron_8>() { return iohelper::HEX1; } +template <> iohelper::ElemType getIOHelperType<_segment_2>() { return iohelper::LINE1; } +template <> iohelper::ElemType getIOHelperType<_segment_3>() { return iohelper::LINE2; } +template <> iohelper::ElemType getIOHelperType<_triangle_3>() { return iohelper::TRIANGLE1; } +template <> iohelper::ElemType getIOHelperType<_triangle_6>() { return iohelper::TRIANGLE2; } +template <> iohelper::ElemType getIOHelperType<_quadrangle_4>() { return iohelper::QUAD1; } +template <> iohelper::ElemType getIOHelperType<_quadrangle_8>() { return iohelper::QUAD2; } +template <> iohelper::ElemType getIOHelperType<_tetrahedron_4>() { return iohelper::TETRA1; } +template <> iohelper::ElemType getIOHelperType<_tetrahedron_10>() { return iohelper::TETRA2; } +template <> iohelper::ElemType getIOHelperType<_hexahedron_8>() { return iohelper::HEX1; } iohelper::ElemType getIOHelperType(ElementType type) { iohelper::ElemType ioh_type = iohelper::MAX_ELEM_TYPE; #define GET_IOHELPER_TYPE(type) \ ioh_type = getIOHelperType(); AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(GET_IOHELPER_TYPE); #undef GET_IOHELPER_TYPE return ioh_type; } /* ------------------------------------------------------------------------ */ void paraviewInit(iohelper::Dumper & dumper, const SolidMechanicsModel & model, const ElementType & type, const std::string & filename) { const Mesh & mesh = model.getFEM().getMesh(); UInt spatial_dimension = mesh.getSpatialDimension(type); UInt nb_nodes = mesh.getNbNodes(); UInt nb_element = mesh.getNbElement(type); std::stringstream filename_sstr; filename_sstr << filename << "_" << type; UInt whoami = StaticCommunicator::getStaticCommunicator()->whoAmI(); UInt nproc = StaticCommunicator::getStaticCommunicator()->getNbProc(); dumper.SetMode(iohelper::TEXT); dumper.SetParallelContext(whoami, nproc); dumper.SetPoints(mesh.getNodes().values, spatial_dimension, nb_nodes, filename_sstr.str().c_str()); dumper.SetConnectivity((int *)mesh.getConnectivity(type).values, getIOHelperType(type), nb_element, iohelper::C_MODE); dumper.AddNodeDataField(model.getDisplacement().values, spatial_dimension, "displacements"); dumper.AddNodeDataField(model.getVelocity().values, spatial_dimension, "velocity"); dumper.AddNodeDataField(model.getAcceleration().values, spatial_dimension, "acceleration"); dumper.AddNodeDataField(model.getMass().values, spatial_dimension, "mass"); dumper.AddNodeDataField(model.getResidual().values, spatial_dimension, "force"); dumper.AddNodeDataField(model.getForce().values, spatial_dimension, "applied_force"); dumper.AddElemDataField(model.getMaterial(0).getStrain(type).values, spatial_dimension*spatial_dimension, "strain"); dumper.AddElemDataField(model.getMaterial(0).getStress(type).values, spatial_dimension*spatial_dimension, "stress"); if(dynamic_cast(&model.getMaterial(0))) { const MaterialDamage & mat = dynamic_cast(model.getMaterial(0)); Real * dam = mat.getDamage(type).storage(); dumper.AddElemDataField(dam, 1, "damage"); } dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("applied_force", 1); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); } /* ------------------------------------------------------------------------ */ void paraviewDump(iohelper::Dumper & dumper) { dumper.Dump(); } /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ // Vector checkpoint_displacements(0, spatial_dimension); // Vector checkpoint_velocity(0, spatial_dimension); // Vector checkpoint_acceleration(0, spatial_dimension); // Vector checkpoint_force(0, spatial_dimension); // /* ------------------------------------------------------------------------ */ // void checkpointInit(iohelper::Dumper & dumper, // const SolidMechanicsModel & model, // const ElementType & type, // const std::string & filename) { // UInt whoami = StaticCommunicator::getStaticCommunicator()->whoAmI(); // UInt nproc = StaticCommunicator::getStaticCommunicator()->getNbProc(); // Vector & displacements = model.getDisplacement(); // Vector & velocity = model.getVelocity(); // Vector & acceleration = model.getAcceleration(); // Vector & force = model.getForce(); // DOFSynchronizer & dof_synchronizer = const_cast(model.getDOFSynchronizer()); // dof_synchronizer.initScatterGatherCommunicationScheme(); // if(whoami == 0){ // const Mesh & mesh = model.getFEM().getMesh(); // UInt nb_nodes = mesh.getNbGlobalNodes(); // checkpoint_displacements.resize(nb_nodes); // checkpoint_velocity .resize(nb_nodes); // checkpoint_acceleration .resize(nb_nodes); // checkpoint_force .resize(nb_nodes); // dof_synchronizer.gather(displacements, 0, &checkpoint_displacements); // dof_synchronizer.gather(velocity , 0, &checkpoint_velocity ); // dof_synchronizer.gather(acceleration , 0, &checkpoint_acceleration ); // dof_synchronizer.gather(force , 0, &checkpoint_force ); // UInt nb_element = mesh.getNbElement(type); // UInt spatial_dimension = mesh.getSpatialDimension(type); // std::stringstream filename_sstr; filename_sstr << filename << "_" << type; // dumper.SetMode(iohelper::COMPRESSED); // dumper.SetParallelContext(whoami, nproc); // dumper.SetPoints(mesh.getNodes().values, // spatial_dimension, nb_nodes, filename_sstr.str().c_str()); // dumper.SetConnectivity((int *)mesh.getConnectivity(type).values, // getIOHelperType(type), nb_element, iohelper::C_MODE); // dumper.AddNodeDataField(checkpoint_displacements.storage(), spatial_dimension, "displacements"); // dumper.AddNodeDataField(checkpoint_velocity .storage(), spatial_dimension, "velocity"); // dumper.AddNodeDataField(checkpoint_acceleration .storage(), spatial_dimension, "acceleration"); // dumper.AddNodeDataField(checkpoint_force .storage(), spatial_dimension, "applied_force"); // dumper.SetPrefix("restart/"); // dumper.Init(); // dumper.Dump(); // } else { // dof_synchronizer.gather(displacements, 0); // dof_synchronizer.gather(velocity , 0); // dof_synchronizer.gather(acceleration , 0); // dof_synchronizer.gather(force , 0); // } // } // /* ------------------------------------------------------------------------ */ // void checkpoint(iohelper::Dumper & dumper, // const SolidMechanicsModel & model) { // UInt whoami = StaticCommunicator::getStaticCommunicator()->whoAmI(); // DOFSynchronizer & dof_synchronizer = const_cast(model.getDOFSynchronizer()); // Vector & displacements = model.getDisplacement(); // Vector & velocity = model.getVelocity(); // Vector & acceleration = model.getAcceleration(); // Vector & force = model.getForce(); // if(whoami == 0){ // dof_synchronizer.gather(displacements, 0, &checkpoint_displacements); // dof_synchronizer.gather(velocity , 0, &checkpoint_velocity ); // dof_synchronizer.gather(acceleration , 0, &checkpoint_acceleration ); // dof_synchronizer.gather(force , 0, &checkpoint_force ); // dumper.Dump(); // } else { // dof_synchronizer.gather(displacements, 0); // dof_synchronizer.gather(velocity , 0); // dof_synchronizer.gather(acceleration , 0); // dof_synchronizer.gather(force , 0); // } // } // /* ------------------------------------------------------------------------ */ // /* ------------------------------------------------------------------------ */ // void restart(const SolidMechanicsModel & model, // const ElementType & type, // const std::string & filename) { // UInt whoami = StaticCommunicator::getStaticCommunicator()->whoAmI(); // UInt nproc = StaticCommunicator::getStaticCommunicator()->getNbProc(); // DOFSynchronizer & dof_synchronizer = const_cast(model.getDOFSynchronizer()); // dof_synchronizer.initScatterGatherCommunicationScheme(); // Vector & displacements = model.getDisplacement(); // Vector & velocity = model.getVelocity(); // Vector & acceleration = model.getAcceleration(); // Vector & force = model.getForce(); // if(whoami == 0){ // const Mesh & mesh = model.getFEM().getMesh(); // UInt nb_nodes = mesh.getNbGlobalNodes(); // UInt spatial_dimension = mesh.getSpatialDimension(type); // std::stringstream filename_sstr; filename_sstr << filename << "_" << type; // iohelper::ReaderRestart reader; // reader.SetMode(iohelper::COMPRESSED); // reader.SetParallelContext(whoami, nproc); // reader.SetPoints(filename_sstr.str().c_str()); // reader.SetConnectivity(getIOHelperType(type)); // reader.AddNodeDataField("displacements"); // reader.AddNodeDataField("velocity"); // reader.AddNodeDataField("acceleration"); // reader.AddNodeDataField("applied_force"); // reader.SetPrefix("restart/"); // reader.Init(); // reader.Read(); // Vector restart_tmp_vect(nb_nodes, spatial_dimension); // displacements.clear(); // velocity.clear(); // acceleration.clear(); // force.clear(); // Real * tmp = reader.GetNodeDataField("displacements"); // std::copy(tmp, tmp + nb_nodes * spatial_dimension, restart_tmp_vect.storage()); // dof_synchronizer.scatter(displacements, 0, &restart_tmp_vect); // tmp = reader.GetNodeDataField("velocity"); // std::copy(tmp, tmp + nb_nodes * spatial_dimension, restart_tmp_vect.storage()); // dof_synchronizer.scatter(velocity , 0, &restart_tmp_vect); // tmp = reader.GetNodeDataField("acceleration"); // std::copy(tmp, tmp + nb_nodes * spatial_dimension, restart_tmp_vect.storage()); // dof_synchronizer.scatter(acceleration , 0, &restart_tmp_vect); // tmp = reader.GetNodeDataField("applied_force"); // std::copy(tmp, tmp + nb_nodes * spatial_dimension, restart_tmp_vect.storage()); // dof_synchronizer.scatter(force , 0, &restart_tmp_vect); // } else { // displacements.clear(); // velocity.clear(); // acceleration.clear(); // force.clear(); // dof_synchronizer.scatter(displacements, 0); // dof_synchronizer.scatter(velocity , 0); // dof_synchronizer.scatter(acceleration , 0); // dof_synchronizer.scatter(force , 0); // } // } /* ------------------------------------------------------------------------ */ diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_non_local.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_material_non_local.cc index bfd055450..dd59d9003 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/test_material_non_local.cc +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_non_local.cc @@ -1,185 +1,185 @@ /** * @file test_non_local_material.cc * @author Nicolas Richart * @date Mon Aug 22 14:07:08 2011 * * @brief test of the main part of the non local materials * * @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 . * */ #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "solid_mechanics_model.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; #ifdef AKANTU_USE_IOHELPER # include "io_helper.hh" static void paraviewInit(iohelper::Dumper & dumper, const SolidMechanicsModel & model); //static void paraviewDump(iohelper::Dumper & dumper); #endif ByElementTypeReal quadrature_points_volumes("quadrature_points_volumes", "test"); const ElementType TYPE = _triangle_6; int main(int argc, char *argv[]) { akantu::initialize(argc, argv); debug::setDebugLevel(akantu::dblWarning); UInt spatial_dimension = 2; Mesh mesh(spatial_dimension); MeshIOMSH mesh_io; mesh_io.read("mesh.msh", mesh); SolidMechanicsModel model(mesh); model.initModel(); model.initVectors(); model.readMaterials("material_non_local.dat"); model.initMaterials(); // model.getFEM().getMesh().initByElementTypeVector(quadrature_points_volumes, 1, 0); const MaterialNonLocal & mat = dynamic_cast &>(model.getMaterial(0)); // mat.computeQuadraturePointsNeighborhoudVolumes(quadrature_points_volumes); Real radius = mat.getRadius(); UInt nb_element = mesh.getNbElement(TYPE); UInt nb_tot_quad = model.getFEM().getNbQuadraturePoints(TYPE) * nb_element; Vector quads(0, spatial_dimension); quads.resize(nb_tot_quad); model.getFEM().interpolateOnQuadraturePoints(mesh.getNodes(), quads, spatial_dimension, TYPE); Vector::iterator first_quad_1 = quads.begin(spatial_dimension); Vector::iterator last_quad_1 = quads.end(spatial_dimension); std::ofstream pout; pout.open("bf_pairs"); UInt q1 = 0; Real R = mat.getRadius(); for(;first_quad_1 != last_quad_1; ++first_quad_1, ++q1) { Vector::iterator first_quad_2 = quads.begin(spatial_dimension); //Vector::iterator last_quad_2 = quads.end(spatial_dimension); UInt q2 = 0; for(;first_quad_2 != last_quad_1; ++first_quad_2, ++q2) { Real d = first_quad_2->distance(*first_quad_1); if(d <= radius) { Real alpha = (1 - d*d/(R*R)); alpha = alpha*alpha; pout << q1 << " " << q2 << " " << alpha << std::endl; } } } pout.close(); mat.savePairs("cl_pairs"); ByElementTypeReal constant("constant_value", "test"); mesh.initByElementTypeVector(constant, 1, 0); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension); for(; it != last_type; ++it) { UInt nb_quadrature_points = model.getFEM().getNbQuadraturePoints(*it); UInt _nb_element = mesh.getNbElement(*it); Vector & constant_vect = constant(*it); constant_vect.resize(_nb_element * nb_quadrature_points); std::fill_n(constant_vect.storage(), nb_quadrature_points * _nb_element, 1.); } ByElementTypeReal constant_avg("constant_value_avg", "test"); mesh.initByElementTypeVector(constant_avg, 1, 0); mat.weightedAvergageOnNeighbours(constant, constant_avg, 1); debug::setDebugLevel(akantu::dblTest); std::cout << constant(TYPE) << std::endl; std::cout << constant_avg(TYPE) << std::endl; debug::setDebugLevel(akantu::dblWarning); #ifdef AKANTU_USE_IOHELPER iohelper::DumperParaview dumper; paraviewInit(dumper, model); #endif akantu::finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ /* iohelper::Dumper vars */ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ template static iohelper::ElemType paraviewType(); -template <> static iohelper::ElemType paraviewType<_segment_2>() { return iohelper::LINE1; } -template <> static iohelper::ElemType paraviewType<_segment_3>() { return iohelper::LINE2; } -template <> static iohelper::ElemType paraviewType<_triangle_3>() { return iohelper::TRIANGLE1; } -template <> static iohelper::ElemType paraviewType<_triangle_6>() { return iohelper::TRIANGLE2; } -template <> static iohelper::ElemType paraviewType<_quadrangle_4>() { return iohelper::QUAD1; } -template <> static iohelper::ElemType paraviewType<_tetrahedron_4>() { return iohelper::TETRA1; } -template <> static iohelper::ElemType paraviewType<_tetrahedron_10>() { return iohelper::TETRA2; } -template <> static iohelper::ElemType paraviewType<_hexahedron_8>() { return iohelper::HEX1; } +template <> iohelper::ElemType paraviewType<_segment_2>() { return iohelper::LINE1; } +template <> iohelper::ElemType paraviewType<_segment_3>() { return iohelper::LINE2; } +template <> iohelper::ElemType paraviewType<_triangle_3>() { return iohelper::TRIANGLE1; } +template <> iohelper::ElemType paraviewType<_triangle_6>() { return iohelper::TRIANGLE2; } +template <> iohelper::ElemType paraviewType<_quadrangle_4>() { return iohelper::QUAD1; } +template <> iohelper::ElemType paraviewType<_tetrahedron_4>() { return iohelper::TETRA1; } +template <> iohelper::ElemType paraviewType<_tetrahedron_10>() { return iohelper::TETRA2; } +template <> iohelper::ElemType paraviewType<_hexahedron_8>() { return iohelper::HEX1; } /* -------------------------------------------------------------------------- */ void paraviewInit(iohelper::Dumper & dumper, const SolidMechanicsModel & model) { UInt spatial_dimension = ElementClass::getSpatialDimension(); UInt nb_nodes = model.getFEM().getMesh().getNbNodes(); UInt nb_element = model.getFEM().getMesh().getNbElement(TYPE); std::stringstream filename; filename << "material_non_local_" << TYPE; dumper.SetMode(iohelper::TEXT); dumper.SetParallelContext(StaticCommunicator::getStaticCommunicator()->whoAmI(), StaticCommunicator::getStaticCommunicator()->getNbProc()); dumper.SetPoints(model.getFEM().getMesh().getNodes().values, spatial_dimension, nb_nodes, filename.str().c_str()); dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(TYPE).values, paraviewType(), nb_element, iohelper::C_MODE); dumper.AddElemDataField(quadrature_points_volumes(TYPE).storage(), 1, "volume"); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); } /* -------------------------------------------------------------------------- */ // void paraviewDump(iohelper::Dumper & dumper) { // dumper.Dump(); // } /* -------------------------------------------------------------------------- */ #endif diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_implicit_dynamic_2d.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_implicit_dynamic_2d.cc index 92e740d02..533eb4eeb 100644 --- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_implicit_dynamic_2d.cc +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_implicit_dynamic_2d.cc @@ -1,315 +1,315 @@ /** * @file test_solid_mechanics_model_implicit_dynamic_2d.cc * @author Nicolas Richart * @date Fri Apr 29 11:32:25 2011 * * @brief test of the dynamic implicit code * * @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 . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "static_communicator.hh" #include "distributed_synchronizer.hh" #include "mesh_partition_scotch.hh" #ifdef AKANTU_USE_SCOTCH #include "mesh_partition_scotch.hh" #endif using namespace akantu; /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER # include "io_helper.hh" static void paraviewInit(iohelper::Dumper & dumper, const SolidMechanicsModel & model); static void paraviewDump(iohelper::Dumper & dumper); #endif /* -------------------------------------------------------------------------- */ const Real F = 0.5e4; #define bar_length 10. #define bar_height 1. #define bar_depth 1. const ElementType TYPE = _triangle_3; UInt spatial_dimension = 2; Real time_step = 1e-4; static Real analytical_solution(Real time) { return 1./pow(M_PI, 4) * ((1. - cos(M_PI*M_PI*time)) + (1. - cos(3*3*M_PI*M_PI*time))/81. + (1. - cos(5*5*M_PI*M_PI*time))/625.); } /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { debug::setDebugLevel(dblWarning); initialize(argc, argv); Mesh mesh(spatial_dimension); StaticCommunicator * comm = StaticCommunicator::getStaticCommunicator(); Int psize = comm->getNbProc(); Int prank = comm->whoAmI(); MeshPartition * partition = NULL; if(prank == 0) { MeshIOMSH mesh_io; mesh_io.read("beam_2d_lin.msh", mesh); partition = new MeshPartitionScotch(mesh, spatial_dimension); partition->reorder(); partition->partitionate(psize); } SolidMechanicsModel * model = new SolidMechanicsModel(mesh); model->initParallel(partition); // UInt nb_nodes = model->getFEM().getMesh().getNbNodes(); /// model initialization model->initVectors(); model->initModel(); model->readMaterials("material_implicit_dynamic.dat"); model->initMaterials(); model->initImplicit(true); // boundary conditions const Vector & position = mesh.getNodes(); Vector & boundary = model->getBoundary(); Vector & force = model->getForce(); Vector & displacment = model->getDisplacement(); //initial conditions model->getForce().clear(); model->getDisplacement().clear(); model->getVelocity().clear(); model->getAcceleration().clear(); // MeshUtils::buildFacets(mesh); // MeshUtils::buildSurfaceID(mesh); // CSR surface_nodes; // MeshUtils::buildNodesPerSurface(mesh, surface_nodes); UInt node_to_print = UInt(-1); bool print_node = false; // for (UInt s = 0; s < surface_nodes.getNbRows(); ++s) { // CSR::iterator snode = surface_nodes.begin(s); // for(; snode != surface_nodes.end(s); ++snode) { // UInt n = *snode; Vector node_to_displace; for (UInt n = 0; n < mesh.getNbNodes(); ++n) { Real x = position(n, 0); Real y = position(n, 1); Real z = 0; if(spatial_dimension == 3) z = position(n, 2); if(Math::are_float_equal(x, 0.) && Math::are_float_equal(y, bar_height / 2.)) { boundary(n,0) = true; boundary(n,1) = true; if(spatial_dimension == 3 && Math::are_float_equal(z, bar_depth / 2.)) boundary(n,2) = true; } if(Math::are_float_equal(x, bar_length) && Math::are_float_equal(y, bar_height / 2.)) { boundary(n,1) = true; if(spatial_dimension == 3 && Math::are_float_equal(z, bar_depth / 2.)) boundary(n,2) = true; } if(Math::are_float_equal(x, bar_length / 2.) && Math::are_float_equal(y, bar_height / 2.)) { if(spatial_dimension < 3 || (spatial_dimension == 3 && Math::are_float_equal(z, bar_depth / 2.))) { force(n,1) = F; if(mesh.isLocalOrMasterNode(n)) { print_node = true; node_to_print = n; std::cout << "I, proc " << prank +1 << " handle the print of node " << n << "(" << x << ", "<< y << ", " << z << ")" << std::endl; } } } } // } model->setTimeStep(time_step); model->updateResidual(); std::stringstream out; out << "position-" << TYPE << "_" << std::scientific << time_step << ".csv"; #ifdef AKANTU_USE_IOHELPER iohelper::DumperParaview dumper; paraviewInit(dumper, *model); #endif std::ofstream pos; if(print_node) { pos.open(out.str().c_str()); if(!pos.good()) { std::cerr << "Cannot open file " << out.str() <assembleStiffnessMatrix(); model->assembleMass(); // model->assembleMassLumped(); // Vector lumped_mass(0,spatial_dimension); // model->getMassMatrix().lump(lumped_mass); // debug::setDebugLevel(dblTest); // std::cout << model->getMass() << lumped_mass; // debug::setDebugLevel(dblWarning); model->getMassMatrix().saveMatrix("M.mtx"); model->getStiffnessMatrix().saveMatrix("K.mtx"); /// time loop for (UInt s = 1; time < 0.62; ++s) { model->implicitPred(); /// convergence loop do { if(count > 0 && prank == 0) std::cout << "passing step " << s << " " << s * time_step << "s - " << std::setw(4) << count << " : " << std::scientific << error << "\r" << std::flush; model->updateResidual(); model->solveDynamic(); model->implicitCorr(); count++; } while(!model->testConvergenceIncrement(1e-12, error) && (count < 1000)); if(prank == 0) std::cout << "passing step " << s << " " << s * time_step << "s - " << std::setw(4) << count << " : " << std::scientific << error << std::endl; count = 0; // if(s % print_freq == 0) { // std::cout << "passing step " << s << "/" << max_steps << " " << s * time_step << "s - " << count / print_freq << std::endl; // count = 0; // } if(print_node) pos << s << "," << s * time_step << "," << displacment(node_to_print, 1) << "," << analytical_solution(s*time_step) << std::endl; #ifdef AKANTU_USE_IOHELPER if(s % 10 == 0) paraviewDump(dumper); #endif time += time_step; } if(print_node) pos.close(); delete model; finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ /* iohelper::Dumper vars */ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ template static iohelper::ElemType paraviewType(); -template <> static iohelper::ElemType paraviewType<_segment_2>() { return iohelper::LINE1; } -template <> static iohelper::ElemType paraviewType<_segment_3>() { return iohelper::LINE2; } -template <> static iohelper::ElemType paraviewType<_triangle_3>() { return iohelper::TRIANGLE1; } -template <> static iohelper::ElemType paraviewType<_triangle_6>() { return iohelper::TRIANGLE2; } -template <> static iohelper::ElemType paraviewType<_quadrangle_4>() { return iohelper::QUAD1; } -template <> static iohelper::ElemType paraviewType<_tetrahedron_4>() { return iohelper::TETRA1; } -template <> static iohelper::ElemType paraviewType<_tetrahedron_10>() { return iohelper::TETRA2; } -template <> static iohelper::ElemType paraviewType<_hexahedron_8>() { return iohelper::HEX1; } +template <> iohelper::ElemType paraviewType<_segment_2>() { return iohelper::LINE1; } +template <> iohelper::ElemType paraviewType<_segment_3>() { return iohelper::LINE2; } +template <> iohelper::ElemType paraviewType<_triangle_3>() { return iohelper::TRIANGLE1; } +template <> iohelper::ElemType paraviewType<_triangle_6>() { return iohelper::TRIANGLE2; } +template <> iohelper::ElemType paraviewType<_quadrangle_4>() { return iohelper::QUAD1; } +template <> iohelper::ElemType paraviewType<_tetrahedron_4>() { return iohelper::TETRA1; } +template <> iohelper::ElemType paraviewType<_tetrahedron_10>() { return iohelper::TETRA2; } +template <> iohelper::ElemType paraviewType<_hexahedron_8>() { return iohelper::HEX1; } /* -------------------------------------------------------------------------- */ void paraviewInit(iohelper::Dumper & dumper, const SolidMechanicsModel & model) { UInt spatial_dimension = ElementClass::getSpatialDimension(); UInt nb_nodes = model.getFEM().getMesh().getNbNodes(); UInt nb_element = model.getFEM().getMesh().getNbElement(TYPE); std::stringstream filename; filename << "dynamic_implicit_beam_" << TYPE; dumper.SetMode(iohelper::TEXT); dumper.SetParallelContext(StaticCommunicator::getStaticCommunicator()->whoAmI(), StaticCommunicator::getStaticCommunicator()->getNbProc()); dumper.SetPoints(model.getFEM().getMesh().getNodes().values, spatial_dimension, nb_nodes, filename.str().c_str()); dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(TYPE).values, paraviewType(), nb_element, iohelper::C_MODE); dumper.AddNodeDataField(model.getDisplacement().values, spatial_dimension, "displacements"); dumper.AddNodeDataField(model.getVelocity().values, spatial_dimension, "velocity"); dumper.AddNodeDataField(model.getAcceleration().values, spatial_dimension, "acceleration"); dumper.AddNodeDataField(model.getResidual().values, spatial_dimension, "residual"); dumper.AddNodeDataField(model.getForce().values, spatial_dimension, "applied_force"); dumper.AddElemDataField(model.getMaterial(0).getStrain(TYPE).values, spatial_dimension*spatial_dimension, "strain"); dumper.AddElemDataField(model.getMaterial(0).getStrain(TYPE).values, spatial_dimension*spatial_dimension, "stress"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("applied_force", 1); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); } /* -------------------------------------------------------------------------- */ void paraviewDump(iohelper::Dumper & dumper) { dumper.Dump(); } /* -------------------------------------------------------------------------- */ #endif diff --git a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2.cc b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2.cc index 1dca28cae..eecad8117 100644 --- a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2.cc +++ b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2.cc @@ -1,231 +1,231 @@ /** * @file test_structural_mechanics_model_bernoulli_beam_2.cc * @author Fabian Barras * @date Thu May 12 16:34:09 2011 * * @brief * * @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 . * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh_struct.hh" #include "structural_mechanics_model.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #define TYPE _bernoulli_beam_2 using namespace akantu; #ifdef AKANTU_USE_IOHELPER # include "io_helper.hh" static void paraviewInit(iohelper::Dumper & dumper, const StructuralMechanicsModel & model); static void paraviewDump(iohelper::Dumper & dumper); #endif //Linear load function static void lin_load(__attribute__ ((unused)) double * position, double * load, __attribute__ ((unused)) Real * normal, __attribute__ ((unused)) UInt surface_id){ memset(load,0,sizeof(Real)*3); load[1]= -60000000; } int main(int argc, char *argv[]){ initialize(argc, argv); Mesh beams(2); debug::setDebugLevel(dblWarning); /* -------------------------------------------------------------------------- */ // Defining the mesh UInt nb_nodes=101; UInt nb_element=nb_nodes-1; UInt nb_nodes_h=101; UInt nb_nodes_v= nb_nodes-nb_nodes_h; Vector & nodes = const_cast &>(beams.getNodes()); nodes.resize(nb_nodes); beams.addConnecticityType(_bernoulli_beam_2); Vector & connectivity = const_cast &>(beams.getConnectivity(_bernoulli_beam_2)); connectivity.resize(nb_element); for(UInt i=0; iinitModel(); StructuralMaterial mat1; mat1.E=2.05e11; mat1.I=0.00128; mat1.A=0.01; model->addMaterial(mat1); model->initVectors(); model->initImplicitSolver(); // const Real M = 3600; // Momentum at 3 Vector & forces = model->getForce(); Vector & displacement = model->getDisplacement(); Vector & boundary = model->getBoundary(); Vector element_material = model->getElementMaterial(_bernoulli_beam_2); forces.clear(); displacement.clear(); model->computeForcesFromFunction(lin_load, akantu::_bft_traction); /**forces(0,2)=-N; forces(nb_nodes-1,2)=N; for (UInt i = 0; i < nb_nodes; ++i) { forces(i,1)=-N/nb_nodes; } */ boundary(0,0) = true; boundary(0,1) = true; //boundary(0,2) = true; //boundary(nb_nodes-1,0) = true; boundary(nb_nodes-1,1) = true; //boundary(nb_nodes-1,2) = true; Real error; model->assembleStiffnessMatrix(); model->getStiffnessMatrix().saveMatrix("Kb.mtx"); UInt count = 0; #ifdef AKANTU_USE_IOHELPER iohelper::DumperParaview dumper; paraviewInit(dumper, *model); #endif do { if(count != 0) std::cerr << count << " - " << error << std::endl; model->updateResidual(); model->solve(); count++; } while (!model->testConvergenceIncrement(1e-10, error) && count < 10); std::cerr << count << " - " << error << std::endl; model->getStiffnessMatrix().saveMatrix("Ka.mtx"); #ifdef AKANTU_USE_IOHELPER paraviewDump(dumper); #endif } /* -------------------------------------------------------------------------- */ /* iohelper::Dumper vars */ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ template static iohelper::ElemType paraviewType(); -template <> static iohelper::ElemType paraviewType<_segment_2>() { return iohelper::LINE1; } -template <> static iohelper::ElemType paraviewType<_segment_3>() { return iohelper::LINE2; } -template <> static iohelper::ElemType paraviewType<_triangle_3>() { return iohelper::TRIANGLE1; } -template <> static iohelper::ElemType paraviewType<_triangle_6>() { return iohelper::TRIANGLE2; } -template <> static iohelper::ElemType paraviewType<_quadrangle_4>() { return iohelper::QUAD1; } -template <> static iohelper::ElemType paraviewType<_tetrahedron_4>() { return iohelper::TETRA1; } -template <> static iohelper::ElemType paraviewType<_tetrahedron_10>() { return iohelper::TETRA2; } -template <> static iohelper::ElemType paraviewType<_hexahedron_8>() { return iohelper::HEX1; } -template <> static iohelper::ElemType paraviewType<_bernoulli_beam_2>(){ return iohelper::LINE1; } +template <> iohelper::ElemType paraviewType<_segment_2>() { return iohelper::LINE1; } +template <> iohelper::ElemType paraviewType<_segment_3>() { return iohelper::LINE2; } +template <> iohelper::ElemType paraviewType<_triangle_3>() { return iohelper::TRIANGLE1; } +template <> iohelper::ElemType paraviewType<_triangle_6>() { return iohelper::TRIANGLE2; } +template <> iohelper::ElemType paraviewType<_quadrangle_4>() { return iohelper::QUAD1; } +template <> iohelper::ElemType paraviewType<_tetrahedron_4>() { return iohelper::TETRA1; } +template <> iohelper::ElemType paraviewType<_tetrahedron_10>() { return iohelper::TETRA2; } +template <> iohelper::ElemType paraviewType<_hexahedron_8>() { return iohelper::HEX1; } +template <> iohelper::ElemType paraviewType<_bernoulli_beam_2>(){ return iohelper::LINE1; } /* -------------------------------------------------------------------------- */ void paraviewInit(iohelper::Dumper & dumper, const StructuralMechanicsModel & model) { UInt spatial_dimension = ElementClass::getSpatialDimension(); UInt nb_nodes = model.getFEM().getMesh().getNbNodes(); UInt nb_element = model.getFEM().getMesh().getNbElement(TYPE); std::stringstream filename; filename << "beam"; dumper.SetMode(iohelper::TEXT); dumper.SetPoints(model.getFEM().getMesh().getNodes().values, spatial_dimension, nb_nodes, filename.str().c_str()); dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(TYPE).values, paraviewType(), nb_element, iohelper::C_MODE); dumper.AddNodeDataField(model.getDisplacement().values, 3, "displacements"); dumper.AddNodeDataField(model.getForce().values, 3, "applied_force"); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); } /* -------------------------------------------------------------------------- */ void paraviewDump(iohelper::Dumper & dumper) { dumper.Dump(); } /* -------------------------------------------------------------------------- */ #endif