diff --git a/packages/damage_non_local.cmake b/packages/damage_non_local.cmake index b8e62e343..62cb182b7 100644 --- a/packages/damage_non_local.cmake +++ b/packages/damage_non_local.cmake @@ -1,84 +1,85 @@ #=============================================================================== # @file 25_damage_non_local.cmake # # @author Nicolas Richart # # @date creation: Fri Jun 15 2012 # @date last modification: Fri Jun 13 2014 # # @brief package description for non-local materials # # @section LICENSE # # Copyright (©) 2014 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 . # #=============================================================================== package_declare(damage_non_local DESCRIPTION "Package for Non-local damage constitutives laws Akantu" DEPENDS lapack) package_declare_sources(damage_non_local model/solid_mechanics/materials/material_damage/material_damage_non_local.hh model/solid_mechanics/materials/material_damage/material_marigo_non_local.hh model/solid_mechanics/materials/material_damage/material_marigo_non_local_inline_impl.cc model/solid_mechanics/materials/material_damage/material_mazars_non_local.cc model/solid_mechanics/materials/material_damage/material_mazars_non_local.hh model/solid_mechanics/materials/material_non_local.hh model/solid_mechanics/materials/material_non_local_includes.hh model/solid_mechanics/materials/material_non_local_inline_impl.cc model/solid_mechanics/materials/weight_functions/base_weight_function.hh model/solid_mechanics/materials/weight_functions/base_weight_function_inline_impl.cc model/solid_mechanics/materials/weight_functions/damaged_weight_function.hh model/solid_mechanics/materials/weight_functions/damaged_weight_function_inline_impl.cc model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function.hh model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function_inline_impl.cc model/solid_mechanics/materials/weight_functions/remove_damaged_with_damage_rate_weight_function.hh model/solid_mechanics/materials/weight_functions/remove_damaged_with_damage_rate_weight_function_inline_impl.cc model/solid_mechanics/materials/weight_functions/stress_based_weight_function.hh model/solid_mechanics/materials/weight_functions/stress_based_weight_function.cc model/solid_mechanics/materials/weight_functions/stress_based_weight_function_inline_impl.cc synchronizer/grid_synchronizer.cc synchronizer/grid_synchronizer.hh model/common/non_local_toolbox/non_local_manager.hh model/common/non_local_toolbox/non_local_manager.cc model/common/non_local_toolbox/non_local_manager_inline_impl.cc model/common/non_local_toolbox/non_local_neighborhood_base.hh model/common/non_local_toolbox/non_local_neighborhood_base.cc model/common/non_local_toolbox/non_local_neighborhood_base_inline_impl.cc model/common/non_local_toolbox/non_local_neighborhood.hh - model/common/non_local_toolbox/non_local_neighborhood.cc + model/common/non_local_toolbox/non_local_neighborhood_tmpl.hh + model/common/non_local_toolbox/non_local_neighborhood_inline_impl.cc ) package_declare_material_infos(damage_non_local LIST AKANTU_DAMAGE_NON_LOCAL_MATERIAL_LIST INCLUDE material_non_local_includes.hh ) package_declare_documentation_files(damage_non_local manual-constitutive-laws-non_local.tex manual-appendix-materials-non-local.tex) package_declare_documentation(damage_non_local "This package activates the non local damage feature of AKANTU" "") diff --git a/src/io/parser/parsable.cc b/src/io/parser/parsable.cc index c8b5c34e4..6fb9b7f2e 100644 --- a/src/io/parser/parsable.cc +++ b/src/io/parser/parsable.cc @@ -1,196 +1,195 @@ /** * @file parsable.cc * * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Tue Sep 02 2014 * * @brief Parsable implementation * * @section LICENSE * * Copyright (©) 2014 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 "parsable.hh" #include "aka_random_generator.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ ParsableParam::ParsableParam(): name(""), description(""), param_type(_pat_internal) { } /* -------------------------------------------------------------------------- */ ParsableParam::ParsableParam(std::string name, std::string description, ParamAccessType param_type) : name(name), description(description), param_type(param_type) { } /* -------------------------------------------------------------------------- */ bool ParsableParam::isWritable() const { return param_type & _pat_writable; } /* -------------------------------------------------------------------------- */ bool ParsableParam::isReadable() const { return param_type & _pat_readable; } /* -------------------------------------------------------------------------- */ bool ParsableParam::isInternal() const { return param_type & _pat_internal; } /* -------------------------------------------------------------------------- */ bool ParsableParam::isParsable() const { return param_type & _pat_parsable; } /* -------------------------------------------------------------------------- */ void ParsableParam::setAccessType(ParamAccessType ptype) { this->param_type = ptype; } /* -------------------------------------------------------------------------- */ void ParsableParam::printself(std::ostream & stream) const { stream << " "; if(isInternal()) stream << "iii"; else { if(isReadable()) stream << "r"; else stream << "-"; if(isWritable()) stream << "w"; else stream << "-"; if(isParsable()) stream << "p"; else stream << "-"; } stream << " "; std::stringstream sstr; sstr << name; UInt width = std::max(int(10 - sstr.str().length()), 0); sstr.width(width); if(description != "") { sstr << " [" << description << "]"; } stream << sstr.str(); width = std::max(int(50 - sstr.str().length()), 0); stream.width(width); stream << " : "; } /* -------------------------------------------------------------------------- */ void ParsableParam::parseParam(const ParserParameter & param) { if(!isParsable()) AKANTU_EXCEPTION("The parameter named " << param.getName() << " is not parsable."); } /* -------------------------------------------------------------------------- */ Parsable::~Parsable() { std::map::iterator it, end; for(it = params.begin(); it != params.end(); ++it){ delete it->second; it->second = NULL; } this->params.clear(); } /* -------------------------------------------------------------------------- */ void Parsable::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); std::map::const_iterator it, end; for(it = params.begin(); it != params.end(); ++it){ stream << space; it->second->printself(stream); } SubSections::const_iterator sit, send; for(sit = sub_sections.begin(); sit != sub_sections.end(); ++sit) { sit->second->printself(stream, indent + 1); } } /* -------------------------------------------------------------------------- */ void Parsable::setParamAccessType(const std::string & name, ParamAccessType ptype) { std::map::iterator it = params.find(name); if(it == params.end()) AKANTU_EXCEPTION("No parameter named " << name << " in parsable."); ParsableParam & param = *(it->second); param.setAccessType(ptype); } /* -------------------------------------------------------------------------- */ void Parsable::registerSubSection(const SectionType & type, const std::string & name, Parsable & sub_section) { SubSectionKey key(type, name); sub_sections[key] = &sub_section; } /* -------------------------------------------------------------------------- */ void Parsable::parseParam(const ParserParameter & in_param) { std::map::iterator it = params.find(in_param.getName()); if(it == params.end()) { if(Parser::parser_permissive) { AKANTU_DEBUG_WARNING("No parameter named " << in_param.getName() << " registered in " << pid << "."); return; } else AKANTU_EXCEPTION("No parameter named " << in_param.getName() << " registered in " << pid << "."); } ParsableParam & param = *(it->second); param.parseParam(in_param); } /* -------------------------------------------------------------------------- */ void Parsable::parseSection(const ParserSection & section) { if(section_type != section.getType()) AKANTU_EXCEPTION("The object " << pid << " is meant to parse section of type " << section_type << ", so it cannot parse section of type " << section.getType()); std::pair params = section.getParameters(); Parser::const_parameter_iterator it = params.first; for (; it != params.second; ++it) { parseParam(*it); } Parser::const_section_iterator sit = section.getSubSections().first; for (; sit != section.getSubSections().second; ++sit) { parseSubSection(*sit); } - } /* -------------------------------------------------------------------------- */ void Parsable::parseSubSection(const ParserSection & section) { SubSectionKey key(section.getType(), section.getName()); SubSections::iterator it = sub_sections.find(key); if(it != sub_sections.end()) { it->second->parseSection(section); } else if(!Parser::parser_permissive) { AKANTU_EXCEPTION("No parsable defined for sub sections of type <" << key.first << "," << key.second << "> in " << pid); } else AKANTU_DEBUG_WARNING("No parsable defined for sub sections of type <" << key.first << "," << key.second << "> in " << pid); } __END_AKANTU__ diff --git a/src/io/parser/parser.hh b/src/io/parser/parser.hh index 95e7c7ed1..efb3ed650 100644 --- a/src/io/parser/parser.hh +++ b/src/io/parser/parser.hh @@ -1,406 +1,408 @@ /** * @file parser.hh * * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Thu Aug 21 2014 * * @brief File parser interface * * @section LICENSE * * Copyright (©) 2014 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 "aka_random_generator.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PARSER_HH__ #define __AKANTU_PARSER_HH__ __BEGIN_AKANTU__ #define AKANTU_SECTION_TYPES \ (global) \ (material) \ (model) \ (mesh) \ (heat) \ (contact) \ (friction) \ - (embedded_interface) \ + (embedded_interface) \ (rules) \ + (neighborhoods) \ (non_local) \ + (weight_function) \ (user) \ (not_defined) #define AKANTU_SECTION_TYPES_PREFIX(elem) BOOST_PP_CAT(_st_, elem) #define AKANTU_SECT_PREFIX(s, data, elem) AKANTU_SECTION_TYPES_PREFIX(elem) enum SectionType { BOOST_PP_SEQ_ENUM(BOOST_PP_SEQ_TRANSFORM(AKANTU_SECT_PREFIX, _, AKANTU_SECTION_TYPES)) }; #undef AKANTU_SECT_PREFIX #define AKANTU_SECTION_TYPE_PRINT_CASE(r, data, elem) \ case AKANTU_SECTION_TYPES_PREFIX(elem): { \ stream << BOOST_PP_STRINGIZE(AKANTU_SECTION_TYPES_PREFIX(elem)); \ break; \ } inline std::ostream & operator <<(std::ostream & stream, SectionType type) { switch(type) { BOOST_PP_SEQ_FOR_EACH(AKANTU_SECTION_TYPE_PRINT_CASE, _, AKANTU_SECTION_TYPES) default: stream << "not a SectionType"; break; } return stream; } #undef AKANTU_SECTION_TYPE_PRINT_CASE enum ParserParameterSearchCxt { _ppsc_current_scope = 0x1, _ppsc_parent_scope = 0x2, _ppsc_current_and_parent_scope = 0x3 }; /* ------------------------------------------------------------------------ */ /* Parameters Class */ /* ------------------------------------------------------------------------ */ class ParserSection; class ParserParameter { public: ParserParameter() : parent_section(NULL), name(std::string()), value(std::string()), dbg_filename(std::string()) {} ParserParameter(const std::string & name, const std::string & value, const ParserSection & parent_section) : parent_section(&parent_section), name(name), value(value), dbg_filename(std::string()) {} ParserParameter(const ParserParameter & param) : parent_section(param.parent_section), name(param.name), value(param.value), dbg_filename(param.dbg_filename), dbg_line(param.dbg_line), dbg_column(param.dbg_column) {} virtual ~ParserParameter() {} const std::string & getName() const { return name; } const std::string & getValue() const { return value; } void setDebugInfo(const std::string & filename, UInt line, UInt column) { dbg_filename = filename; dbg_line = line; dbg_column = column; } template inline operator T() const; // template inline operator Vector() const; // template inline operator Matrix() const; void printself(std::ostream & stream, unsigned int indent = 0) const { stream << name << ": " << value << " (" << dbg_filename << ":" << dbg_line << ":" << dbg_column << ")"; } private: void setParent(const ParserSection & sect) { parent_section = § } friend class ParserSection; private: const ParserSection * parent_section; std::string name; std::string value; std::string dbg_filename; UInt dbg_line, dbg_column; }; /* ------------------------------------------------------------------------ */ /* Sections Class */ /* ------------------------------------------------------------------------ */ class ParserSection { public: typedef std::multimap SubSections; typedef std::map Parameters; private: typedef SubSections::const_iterator const_section_iterator_; public: /* ------------------------------------------------------------------------ */ /* SubSection iterator */ /* ------------------------------------------------------------------------ */ class const_section_iterator { public: const_section_iterator(const const_section_iterator & other) : it(other.it) { } const_section_iterator(const const_section_iterator_ & it) : it(it) { } const_section_iterator & operator=(const const_section_iterator & other) { if(this != &other) { it = other.it; } return *this; } const ParserSection & operator* () const { return it->second; } const ParserSection * operator-> () const { return &(it->second); } bool operator==(const const_section_iterator & other) const { return it == other.it; } bool operator!=(const const_section_iterator & other) const { return it != other.it; } const_section_iterator & operator++() { ++it; return *this; } const_section_iterator operator++(int) { const_section_iterator tmp = *this; operator++(); return tmp; } private: const_section_iterator_ it; }; /* ------------------------------------------------------------------------ */ /* Parameters iterator */ /* ------------------------------------------------------------------------ */ class const_parameter_iterator { public: const_parameter_iterator(const const_parameter_iterator & other) : it(other.it) { } const_parameter_iterator(const Parameters::const_iterator & it) : it(it) { } const_parameter_iterator & operator=(const const_parameter_iterator & other) { if(this != &other) { it = other.it; } return *this; } const ParserParameter & operator* () const { return it->second; } const ParserParameter * operator->() { return &(it->second); }; bool operator==(const const_parameter_iterator & other) const { return it == other.it; } bool operator!=(const const_parameter_iterator & other) const { return it != other.it; } const_parameter_iterator & operator++() { ++it; return *this; } const_parameter_iterator operator++(int) { const_parameter_iterator tmp = *this; operator++(); return tmp; } private: Parameters::const_iterator it; }; /* ---------------------------------------------------------------------- */ ParserSection() : parent_section(NULL), name(std::string()), type(_st_not_defined){} ParserSection(const std::string & name, SectionType type) : parent_section(NULL), name(name), type(type) {} ParserSection(const std::string & name, SectionType type, const std::string & option, const ParserSection & parent_section) : parent_section(&parent_section), name(name), type(type), option(option) {} ParserSection(const ParserSection & section) : parent_section(section.parent_section), name(section.name), type(section.type), option(section.option), parameters(section.parameters), sub_sections_by_type(section.sub_sections_by_type) { setChldrenPointers(); } ParserSection & operator=(const ParserSection & other) { if(&other != this) { parent_section = other.parent_section; name = other.name; type = other.type; option = other.option; parameters = other.parameters; - sub_sections_by_type = sub_sections_by_type; + sub_sections_by_type = other.sub_sections_by_type; setChldrenPointers(); } return *this; } virtual ~ParserSection(); virtual void printself(std::ostream & stream, unsigned int indent = 0) const; /* ---------------------------------------------------------------------- */ /* Creation functions */ /* ---------------------------------------------------------------------- */ public: ParserParameter & addParameter(const ParserParameter & param); ParserSection & addSubSection(const ParserSection & section); private: void setChldrenPointers() { Parameters::iterator pit = this->parameters.begin(); for(;pit != this->parameters.end(); ++pit) pit->second.setParent(*this); SubSections::iterator sit = this->sub_sections_by_type.begin(); for (;sit != this->sub_sections_by_type.end(); ++sit) sit->second.setParent(*this); } /* ---------------------------------------------------------------------- */ /* Accessors */ /* ---------------------------------------------------------------------- */ public: std::pair getSubSections(SectionType type = _st_not_defined) const { if(type != _st_not_defined) { std::pair range = sub_sections_by_type.equal_range(type); return std::pair(range.first, range.second); } else { return std::pair(sub_sections_by_type.begin(), sub_sections_by_type.end()); } } std::pair getParameters() const { return std::pair(parameters.begin(), parameters.end()); } /* ---------------------------------------------------------------------- */ const ParserParameter & getParameter(const std::string & name, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { Parameters::const_iterator it; if(search_ctx & _ppsc_current_scope) it = parameters.find(name); if(it == parameters.end()) { if((search_ctx & _ppsc_parent_scope) && parent_section) return parent_section->getParameter(name, search_ctx); else { AKANTU_SILENT_EXCEPTION("The parameter " << name << " has not been found in the specified context"); } } return it->second; } /* ------------------------------------------------------------------------ */ bool hasParameter(const std::string & name, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { Parameters::const_iterator it; if(search_ctx & _ppsc_current_scope) it = parameters.find(name); if(it == parameters.end()) { if((search_ctx & _ppsc_parent_scope) && parent_section) return parent_section->hasParameter(name, search_ctx); else { return false; } } return true; } /* -------------------------------------------------------------------------- */ template T getParameterValue(const std::string & name, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { const ParserParameter & tmp_param = getParameter(name, search_ctx); T t = tmp_param; return t; } /* -------------------------------------------------------------------------- */ const std::string & getName() const { return name; } const SectionType & getType() const { return type; } const std::string & getOption() const { return option; } protected: void setParent(const ParserSection & sect) { parent_section = § } /* ---------------------------------------------------------------------- */ /* Members */ /* ---------------------------------------------------------------------- */ private: const ParserSection * parent_section; std::string name; SectionType type; std::string option; Parameters parameters; SubSections sub_sections_by_type; }; /* ------------------------------------------------------------------------ */ /* Parser Class */ /* ------------------------------------------------------------------------ */ class Parser : public ParserSection { public: Parser() : ParserSection("global", _st_global) {} void parse(const std::string & filename); std::string getLastParsedFile() const; static bool isPermissive() { return parser_permissive; } public: static Real parseReal (const std::string & value, const ParserSection & section); static Vector parseVector(const std::string & value, const ParserSection & section); static Matrix parseMatrix(const std::string & value, const ParserSection & section); static RandomParameter parseRandomParameter(const std::string & value, const ParserSection & section); protected: template static T parseType(const std::string & value, Grammar & grammar); protected: friend class Parsable; static bool parser_permissive; std::string last_parsed_file; }; inline std::ostream & operator <<(std::ostream & stream, const ParserParameter &_this) { _this.printself(stream); return stream; } inline std::ostream & operator <<(std::ostream & stream, const ParserSection & section) { section.printself(stream); return stream; } __END_AKANTU__ #include "parser_tmpl.hh" #endif /* __AKANTU_PARSER_HH__ */ diff --git a/src/model/common/non_local_toolbox/non_local_manager.cc b/src/model/common/non_local_toolbox/non_local_manager.cc index af43d97fe..b9c701702 100644 --- a/src/model/common/non_local_toolbox/non_local_manager.cc +++ b/src/model/common/non_local_toolbox/non_local_manager.cc @@ -1,181 +1,295 @@ /** * @file non_local_manager.cc * @author Aurelia Isabel Cuba Ramos * @date Mon Sep 21 15:32:10 2015 * * @brief Implementation of non-local manager * * @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 "non_local_manager.hh" #include "non_local_neighborhood.hh" #include "material_non_local.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ NonLocalManager::NonLocalManager(SolidMechanicsModel & model, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), + Parsable(_st_neighborhoods, id), model(model), default_neighborhood("default_neighborhood"), quad_positions("quad_positions", id), - volumes("volumes", id) { - UInt spatial_dimension = this->model.getSpatialDimension(); + volumes("volumes", id), + spatial_dimension(this->model.getSpatialDimension()), + compute_stress_calls(0){ Mesh & mesh = this->model.getMesh(); /// initialize the element type map arrays - mesh.initElementTypeMapArray(volumes, 1, spatial_dimension, false, _ek_regular, true); mesh.initElementTypeMapArray(quad_positions, spatial_dimension, spatial_dimension, false, _ek_regular, true); -#ifdef AKANTU_IGFEM - mesh.initElementTypeMapArray(volumes, 1, spatial_dimension, false, _ek_igfem, true); - mesh.initElementTypeMapArray(quad_positions, spatial_dimension, spatial_dimension, false, _ek_igfem, true); -#endif + /// parse the neighborhood information from the input file + const Parser & parser = getStaticParser(); + const ParserSection & non_local_section = *(parser.getSubSections(_st_neighborhoods).first); + /// iterate over all the non-local sections and store them in a map + std::pair + weight_sect = non_local_section.getSubSections(_st_non_local); + Parser::const_section_iterator it = weight_sect.first; + for (; it != weight_sect.second; ++it) { + const ParserSection & section = *it; + ID name = section.getName(); + this->weight_function_types[name] = section; + } } /* -------------------------------------------------------------------------- */ NonLocalManager::~NonLocalManager() { NeighborhoodMap::iterator it; for (it = neighborhoods.begin(); it != neighborhoods.end(); ++it) { if(it->second) delete it->second; } } /* -------------------------------------------------------------------------- */ void NonLocalManager::setJacobians(const FEEngine & fe_engine, const ElementKind & kind) { Mesh & mesh = this->model.getMesh(); - UInt spatial_dimension = this->model.getSpatialDimension(); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, kind); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, gt, kind); for(; it != last_type; ++it) { jacobians(*it, gt) = &fe_engine.getIntegratorInterface().getJacobians(*it, gt); } } } /* -------------------------------------------------------------------------- */ -void NonLocalManager::createNeighborhood(const ID & type, Real radius, const ID & name) { +void NonLocalManager::createNeighborhood(const ID & weight_func, const ID & neighborhood_id) { AKANTU_DEBUG_IN(); - /// check if neighborhood already exists - NeighborhoodMap::const_iterator it = neighborhoods.find(name); + + const ParserSection & section = this->weight_function_types[weight_func]; + const ID weight_func_type = section.getOption(); + /// create new neighborhood for given ID + std::stringstream sstr; sstr << id << ":neighborhood:" << neighborhood_id; - if (it == neighborhoods.end() ) { - /// create new neighborhood for given ID - std::stringstream sstr; sstr << "neighborhood:" << name; - neighborhoods[name] = new NonLocalNeighborhood(*this, radius, type, this->quad_positions, sstr.str()); - } + if (weight_func_type == "base_wf") + neighborhoods[neighborhood_id] = new NonLocalNeighborhood(*this, this->quad_positions, sstr.str()); + else + {std::cout << "error here" << std::endl;} + + neighborhoods[neighborhood_id]->parseSection(section); + neighborhoods[neighborhood_id]->initNeighborhood(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NonLocalManager::createNeighborhoodSynchronizers() { NeighborhoodMap::iterator it; for (it = neighborhoods.begin(); it != neighborhoods.end(); ++it) { it->second->createGridSynchronizer(); } } /* -------------------------------------------------------------------------- */ void NonLocalManager::flattenInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) { - - internal_flat.clear(); const ID field_name = internal_flat.getName(); for (UInt m = 0; m < this->non_local_materials.size(); ++m) { Material & material = *(this->non_local_materials[m]); - if (material.isInternal(field_name, kind)) + if (material.isInternal(field_name, kind)) material.flattenInternal(field_name, internal_flat, ghost_type, kind); } } /* -------------------------------------------------------------------------- */ void NonLocalManager::averageInternals(const GhostType & ghost_type) { - /// update the flattened version of the internals - std::map::iterator non_local_variable_it = non_local_variables.begin(); - std::map::iterator non_local_variable_end = non_local_variables.end(); - - for(; non_local_variable_it != non_local_variable_end; ++non_local_variable_it) { - for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { - GhostType ghost_type = (GhostType) gt; - this->flattenInternal(non_local_variable_it->second->local, ghost_type, _ek_regular); - this->flattenInternal(non_local_variable_it->second->non_local, ghost_type, _ek_regular); - } - } - /// loop over all neighborhoods and compute the non-local variables NeighborhoodMap::iterator neighborhood_it = neighborhoods.begin(); - NeighborhoodMap::iterator neighborhood_end = neighborhoods.end(); + NeighborhoodMap::iterator neighborhood_end = neighborhoods.end(); for (; neighborhood_it != neighborhood_end; ++neighborhood_it) { - non_local_variable_it = non_local_variables.begin(); - non_local_variable_end = non_local_variables.end(); + /// update the weights of the weight function + if (ghost_type == _not_ghost) + neighborhood_it->second->updateWeights(); + /// loop over all the non-local variables + std::map::iterator non_local_variable_it = non_local_variables.begin(); + std::map::iterator non_local_variable_end = non_local_variables.end(); for(; non_local_variable_it != non_local_variable_end; ++non_local_variable_it) { NonLocalVariable * non_local_var = non_local_variable_it->second; - neighborhood_it->second->weightedAvergageOnNeighbours(non_local_var->local, non_local_var->non_local, non_local_var->nb_component, ghost_type); + neighborhood_it->second->weightedAverageOnNeighbours(non_local_var->local, non_local_var->non_local, non_local_var->nb_component, ghost_type); + } - } + } + } /* -------------------------------------------------------------------------- */ void NonLocalManager::init(){ - UInt spatial_dimension = this->model.getSpatialDimension(); - /// exchange the missing ghosts for the non-local neighborhoods this->createNeighborhoodSynchronizers(); /// insert the ghost quadrature points of the non-local materials into the non-local neighborhoods for(UInt m = 0; m < this->non_local_materials.size(); ++m) { switch (spatial_dimension) { case 1: dynamic_cast &>(*(this->non_local_materials[m])).insertQuadsInNeighborhoods(_ghost); break; case 2: dynamic_cast &>(*(this->non_local_materials[m])).insertQuadsInNeighborhoods(_ghost); break; case 3: dynamic_cast &>(*(this->non_local_materials[m])).insertQuadsInNeighborhoods(_ghost); break; } } this->setJacobians(this->model.getFEEngine(), _ek_regular); + this->initNonLocalVariables(); this->updatePairLists(); + this->initElementTypeMap(1, volumes); this->computeWeights(); } +/* -------------------------------------------------------------------------- */ +void NonLocalManager::initNonLocalVariables(){ + /// loop over all the non-local variables + std::map::iterator non_local_variable_it = non_local_variables.begin(); + std::map::iterator non_local_variable_end = non_local_variables.end(); + + for(; non_local_variable_it != non_local_variable_end; ++non_local_variable_it) { + NonLocalVariable & variable = *(non_local_variable_it->second); + this->initElementTypeMap(variable.nb_component, variable.non_local); + } +} + +/* -------------------------------------------------------------------------- */ +void NonLocalManager::initElementTypeMap(UInt nb_component, ElementTypeMapReal & element_map) { + Mesh & mesh = this->model.getMesh(); + /// need to resize the arrays + for(UInt g = _not_ghost; g <= _ghost; ++g) { + GhostType gt = (GhostType) g; + Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_regular); + Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_regular); + for(; it != end; ++it) { + ElementType el_type = *it; + UInt nb_element = mesh.getNbElement(*it, gt); + UInt nb_quads = this->model.getFEEngine().getNbQuadraturePoints(*it, gt); + if (!element_map.exists(el_type, gt)) { + element_map.alloc(nb_element * nb_quads, + nb_component, + el_type, + gt); + } + } + } +} + +/* -------------------------------------------------------------------------- */ +void NonLocalManager::distributeInternals(ElementKind kind) { + + /// loop over all the non-local variables and copy back their values into the materials + std::map::iterator non_local_variable_it = non_local_variables.begin(); + std::map::iterator non_local_variable_end = non_local_variables.end(); + for(; non_local_variable_it != non_local_variable_end; ++non_local_variable_it) { + NonLocalVariable * non_local_var = non_local_variable_it->second; + const ID field_name = non_local_var->non_local.getName(); + /// loop over all the materials + for (UInt m = 0; m < this->non_local_materials.size(); ++m) { + if (this->non_local_materials[m]->isInternal(field_name, kind)) + + switch (spatial_dimension) { + case 1: + dynamic_cast &>(*(this->non_local_materials[m])).updateNonLocalInternals(non_local_var->non_local, field_name, non_local_var->nb_component); + break; + case 2: + dynamic_cast &>(*(this->non_local_materials[m])).updateNonLocalInternals(non_local_var->non_local, field_name, non_local_var->nb_component); + break; + case 3: + dynamic_cast &>(*(this->non_local_materials[m])).updateNonLocalInternals(non_local_var->non_local, field_name, non_local_var->nb_component); + break; + } + } + } +} + +/* -------------------------------------------------------------------------- */ +void NonLocalManager::computeAllNonLocalStresses() { + + /// update the flattened version of the internals + std::map::iterator non_local_variable_it = non_local_variables.begin(); + std::map::iterator non_local_variable_end = non_local_variables.end(); + + for(; non_local_variable_it != non_local_variable_end; ++non_local_variable_it) { + non_local_variable_it->second->local.clear(); + non_local_variable_it->second->non_local.clear(); + for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { + GhostType ghost_type = (GhostType) gt; + this->flattenInternal(non_local_variable_it->second->local, ghost_type, _ek_regular); + } + } + + this->volumes.clear(); + SynchronizerRegistry & synch_registry = this->model.getSynchronizerRegistry(); + synch_registry.asynchronousSynchronize(_gst_mnl_for_average); + + AKANTU_DEBUG_INFO("Compute non local variables for local elements"); + synch_registry.asynchronousSynchronize(_gst_mnl_weight); + this->averageInternals(_not_ghost); + + AKANTU_DEBUG_INFO("Wait distant non local stresses"); + synch_registry.waitEndSynchronize(_gst_mnl_for_average); + + this->averageInternals(_ghost); + + /// copy the results in the materials + this->distributeInternals(_ek_regular); + /// loop over all the materials and update the weights + for (UInt m = 0; m < this->non_local_materials.size(); ++m) { + switch (spatial_dimension) { + case 1: + dynamic_cast &>(*(this->non_local_materials[m])).computeNonLocalStresses(_not_ghost); break; + case 2: + dynamic_cast &>(*(this->non_local_materials[m])).computeNonLocalStresses(_not_ghost); break; + case 3: + dynamic_cast &>(*(this->non_local_materials[m])).computeNonLocalStresses(_not_ghost); break; + } + } + + ++this->compute_stress_calls; +} + + __END_AKANTU__ diff --git a/src/model/common/non_local_toolbox/non_local_manager.hh b/src/model/common/non_local_toolbox/non_local_manager.hh index bfa10438b..c52b832c1 100644 --- a/src/model/common/non_local_toolbox/non_local_manager.hh +++ b/src/model/common/non_local_toolbox/non_local_manager.hh @@ -1,160 +1,196 @@ /** * @file non_local_manager.hh * @author Aurelia Isabel Cuba Ramos * @date Mon Sep 21 14:21:33 2015 * * @brief Classes that manages all the non-local neighborhoods * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_NON_LOCAL_MANAGER_HH__ #define __AKANTU_NON_LOCAL_MANAGER_HH__ /* -------------------------------------------------------------------------- */ #include "aka_memory.hh" #include "solid_mechanics_model.hh" #include "non_local_neighborhood_base.hh" +#include "parsable.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ -class NonLocalManager : public Memory { +class NonLocalManager : public Memory, public Parsable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NonLocalManager(SolidMechanicsModel & model, const ID & id = "non_local_manager", const MemoryID & memory_id = 0); virtual ~NonLocalManager(); typedef std::map NeighborhoodMap; typedef std::pair KeyCOO; /* -------------------------------------------------------------------------- */ /* Methods */ /* -------------------------------------------------------------------------- */ public: /// initialize the non-local manager: compute pair lists and weights for all neighborhoods void init(); /// insert new quadrature point in the grid - inline void insertQuad(const QuadraturePoint & quad, const Vector & coords, Real radius, const ID & type, ID name = ""); + inline void insertQuad(const QuadraturePoint & quad, const Vector & coords, const ID & weight_func, ID neighborhood = ""); /// return the fem object associated with a provided name inline NonLocalNeighborhoodBase & getNeighborhood(const ID & name = "") const; /// create the grid synchronizers for each neighborhood void createNeighborhoodSynchronizers(); /// compute the weights in each neighborhood for non-local averaging inline void computeWeights(); /// compute the weights in each neighborhood for non-local averaging inline void updatePairLists(); /// register a new non-local material inline void registerNonLocalMaterial(Material & new_mat); /// register a non-local variable inline void registerNonLocalVariable(const ID & variable_name, const ID & nl_variable_name, UInt nb_component); + /// average the non-local variables void averageInternals(const GhostType & ghost_type = _not_ghost); + /// average the non-local variables + //void testo(const GhostType & ghost_type = _not_ghost); + + /// average the internals and compute the non-local stresses + void computeAllNonLocalStresses(); + + /// register a new internal needed for the weight computations + inline ElementTypeMapReal & registerWeightFunctionInternal(const ID & field_name); + + /// update the flattened version of the weight function internals + inline void updateWeightFunctionInternals(); + private: /// create a new neighborhood for a given domain ID - void createNeighborhood(const ID & type, Real radius, const ID & name); + void createNeighborhood(const ID & weight_func, const ID & neighborhood); /// flatten the material internal fields needed for the non-local computations void flattenInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind); /// set the values of the jacobians void setJacobians(const FEEngine & fe_engine, const ElementKind & kind); + /// allocation of elment type maps + void initElementTypeMap(UInt nb_component, ElementTypeMapReal & element_map); + + /// allocate the non-local variables + void initNonLocalVariables(); + + /// copy the results of the averaging in the materials + void distributeInternals(ElementKind kind); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: + AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); AKANTU_GET_MACRO(Model, model, const SolidMechanicsModel &); AKANTU_GET_MACRO_NOT_CONST(Volumes, volumes, ElementTypeMapReal &) + AKANTU_GET_MACRO(NbStressCalls, compute_stress_calls, UInt); inline const Array & getJacobians(const ElementType & type, const GhostType & ghost_type) { return *jacobians(type, ghost_type); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// the non-local neighborhoods present NeighborhoodMap neighborhoods; /// list of all the non-local materials in the model std::vector non_local_materials; struct NonLocalVariable { NonLocalVariable(const ID & variable_name, const ID & nl_variable_name, const ID & id, UInt nb_component) : local(variable_name, id), non_local(nl_variable_name, id), nb_component(nb_component){ } ElementTypeMapReal local; ElementTypeMapReal non_local; UInt nb_component; }; /// the non-local variables associated to a certain neighborhood std::map non_local_variables; /// reference to the model SolidMechanicsModel & model; /// jacobians for all the elements in the mesh ElementTypeMap * > jacobians; /// default neighborhood object std::string default_neighborhood; /// store the position of the quadrature points ElementTypeMapReal quad_positions; /// store the volume of each quadrature point for the non-local weight normalization ElementTypeMapReal volumes; + /// the spatial dimension + const UInt spatial_dimension; + + /// counter for computeStress calls + UInt compute_stress_calls; + + /// map to store weight function types from input file + std::map weight_function_types; + + /// map to store the internals needed by the weight functions + std::map weight_function_internals; + }; __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "non_local_manager_inline_impl.cc" #endif /* __AKANTU_NON_LOCAL_MANAGER_HH__ */ diff --git a/src/model/common/non_local_toolbox/non_local_manager_inline_impl.cc b/src/model/common/non_local_toolbox/non_local_manager_inline_impl.cc index cb734ba1a..68324929a 100644 --- a/src/model/common/non_local_toolbox/non_local_manager_inline_impl.cc +++ b/src/model/common/non_local_toolbox/non_local_manager_inline_impl.cc @@ -1,145 +1,172 @@ /** * @file non_local_manager_inline_impl.cc * @author Aurelia Isabel Cuba Ramos * @date Mon Sep 21 17:30:03 2015 * * @brief inline implementation of non-local manager functions * * @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 "grid_synchronizer.hh" #include "synchronizer_registry.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_NON_LOCAL_MANAGER_INLINE_IMPL_CC__ #define __AKANTU_NON_LOCAL_MANAGER_INLINE_IMPL_CC__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline void NonLocalManager::insertQuad(const QuadraturePoint & quad, const Vector & coords, - Real radius, const ID & type, ID name) { + const ID & weight_func, ID neighborhood) { - if (name == "") name = default_neighborhood; + if (neighborhood == "") neighborhood = default_neighborhood; - NeighborhoodMap::const_iterator it = neighborhoods.find(name); + NeighborhoodMap::const_iterator it = neighborhoods.find(neighborhood); if (it == neighborhoods.end()) { - this->createNeighborhood(type, radius, name); + this->createNeighborhood(weight_func, neighborhood); } - neighborhoods[name]->insertQuad(quad, coords); + neighborhoods[neighborhood]->insertQuad(quad, coords); } /* -------------------------------------------------------------------------- */ inline NonLocalNeighborhoodBase & NonLocalManager::getNeighborhood(const ID & name) const{ AKANTU_DEBUG_IN(); ID tmp_name = name; if (name == "") tmp_name = default_neighborhood; NeighborhoodMap::const_iterator it = neighborhoods.find(tmp_name); AKANTU_DEBUG_ASSERT(it != neighborhoods.end(), "The neighborhood " << tmp_name << " is not registered"); AKANTU_DEBUG_OUT(); return *(it->second); } /* -------------------------------------------------------------------------- */ inline void NonLocalManager::computeWeights() { AKANTU_DEBUG_IN(); - Mesh & mesh = this->model.getMesh(); - UInt spatial_dimension = this->model.getSpatialDimension(); - /// need to resize the arrays - for(UInt g = _not_ghost; g <= _ghost; ++g) { - GhostType gt = (GhostType) g; - Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_regular); - Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_regular); - for(; it != end; ++it) { - UInt nb_element = mesh.getNbElement(*it, gt); - UInt nb_quads = this->model.getFEEngine().getNbQuadraturePoints(*it, gt); - volumes(*it, gt).resize(nb_element *nb_quads); - } - } - + // /// need to resize the arrays + // for(UInt g = _not_ghost; g <= _ghost; ++g) { + // GhostType gt = (GhostType) g; + // Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_regular); + // Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_regular); + // for(; it != end; ++it) { + // UInt nb_element = mesh.getNbElement(*it, gt); + // UInt nb_quads = this->model.getFEEngine().getNbQuadraturePoints(*it, gt); + // volumes(*it, gt).resize(nb_element *nb_quads); + // } + // } + + this->updateWeightFunctionInternals(); this->volumes.clear(); NeighborhoodMap::iterator it = neighborhoods.begin(); NeighborhoodMap::iterator end = neighborhoods.end(); for (; it != end; ++it) it->second->computeWeights(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void NonLocalManager::updatePairLists() { AKANTU_DEBUG_IN(); /// compute the position of the quadrature points this->model.getFEEngine().computeQuadraturePointsCoordinates(quad_positions); NeighborhoodMap::iterator it = neighborhoods.begin(); NeighborhoodMap::iterator end = neighborhoods.end(); for (; it != end; ++it) it->second->updatePairList(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void NonLocalManager::registerNonLocalVariable(const ID & variable_name, const ID & nl_variable_name, UInt nb_component) { AKANTU_DEBUG_IN(); std::map::iterator non_local_variable_it = non_local_variables.find(variable_name); if (non_local_variable_it == non_local_variables.end()) non_local_variables[variable_name] = new NonLocalVariable(variable_name, nl_variable_name, this->id, nb_component); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void NonLocalManager::registerNonLocalMaterial(Material & new_mat) { non_local_materials.push_back(&new_mat); } +/* -------------------------------------------------------------------------- */ +inline ElementTypeMapReal & NonLocalManager::registerWeightFunctionInternal(const ID & field_name) { + + AKANTU_DEBUG_IN(); + + std::map::const_iterator it = this->weight_function_internals.find(field_name); + if (it == weight_function_internals.end()) { + weight_function_internals[field_name] = new ElementTypeMapReal(); + } + + return *(weight_function_internals[field_name]); + AKANTU_DEBUG_OUT(); + +} + +/* -------------------------------------------------------------------------- */ +inline void NonLocalManager::updateWeightFunctionInternals() { + + std::map::const_iterator it = this->weight_function_internals.begin(); + std::map::const_iterator end = this->weight_function_internals.end(); + for (; it != end; ++it) { + it->second->clear(); + for(UInt g = _not_ghost; g <= _ghost; ++g) { + GhostType ghost_type = (GhostType) g; + this->flattenInternal(*(it->second), ghost_type, _ek_regular); + } + } +} __END_AKANTU__ #endif /* __AKANTU_NON_LOCAL_MANAGER_INLINE_IMPL_CC__ */ diff --git a/src/model/common/non_local_toolbox/non_local_neighborhood.hh b/src/model/common/non_local_toolbox/non_local_neighborhood.hh index e4fc30930..dc840a652 100644 --- a/src/model/common/non_local_toolbox/non_local_neighborhood.hh +++ b/src/model/common/non_local_toolbox/non_local_neighborhood.hh @@ -1,128 +1,211 @@ /** * @file non_local_neighborhood.hh * @author Aurelia Isabel Cuba Ramos * @date Sat Sep 26 18:29:59 2015 * * @brief Non-local neighborhood for non-local averaging based on * weight function * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_NON_LOCAL_NEIGHBORHOOD_HH__ #define __AKANTU_NON_LOCAL_NEIGHBORHOOD_HH__ /* -------------------------------------------------------------------------- */ #include "base_weight_function.hh" #include "non_local_neighborhood_base.hh" +#include "parsable.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class NonLocalManager; + class TestWeightFunction; } __BEGIN_AKANTU__ -class WeightFunction { +template +class NonLocalNeighborhood : public NonLocalNeighborhoodBase { + /* ------------------------------------------------------------------------ */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ */ public: + + NonLocalNeighborhood(NonLocalManager & manager, + const ElementTypeMapReal & quad_coordinates, + const ID & id = "neighborhood", + const MemoryID & memory_id = 0); + virtual ~NonLocalNeighborhood(); + + /* ------------------------------------------------------------------------ */ + /* Methods */ + /* ------------------------------------------------------------------------ */ +public: + + /// compute the weights for non-local averaging + void computeWeights(); + + /// save the pair of weights in a file + void saveWeights(const std::string & filename) const; + + /// compute the non-local counter part for a given element type map + virtual void weightedAverageOnNeighbours(const ElementTypeMapReal & to_accumulate, + ElementTypeMapReal & accumulated, + UInt nb_degree_of_freedom, + const GhostType & ghost_type2) const; + + /// update the weights based on the weight function + void updateWeights(); + +protected: + virtual inline UInt getNbDataForElements(const Array & elements, + SynchronizationTag tag) const; + + virtual inline void packElementData(CommunicationBuffer & buffer, + const Array & elements, + SynchronizationTag tag) const; + + virtual inline void unpackElementData(CommunicationBuffer & buffer, + const Array & elements, + SynchronizationTag tag); + + /* -------------------------------------------------------------------------- */ + /* Accessor */ + /* -------------------------------------------------------------------------- */ + AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, const NonLocalManager &); + AKANTU_GET_MACRO_NOT_CONST(NonLocalManager, *non_local_manager, NonLocalManager &); + + /* ------------------------------------------------------------------------ */ + /* Class Members */ + /* ------------------------------------------------------------------------ */ +private: + + /// Pointer to non-local manager class + NonLocalManager * non_local_manager; + + /// the weights associated to the pairs + Array * pair_weight[2]; + + /// weight function + WeightFunction * weight_function; + +}; + + +class TestWeightFunction : public Parsable { + +public: + TestWeightFunction(NonLocalManager & manager) : Parsable(_st_weight_function, "weight_function:"), + manager(manager) +{ + this->registerParam("update_rate" , update_rate, 0U , + _pat_parsmod, "Update frequency"); + + } + void setRadius(Real radius) { /// set the non-local radius and update R^2 accordingly this->R = radius; this->R2 = this->R * this->R; } Real operator()(Real r, const __attribute__((unused)) QuadraturePoint & q1, const __attribute__((unused)) QuadraturePoint & q2) { /// initialize the weight Real w = 0; /// compute weight for given r if(r <= this->R) { Real alpha = (1. - r*r / this->R2); w = alpha * alpha; // *weight = 1 - sqrt(r / radius); } return w; } void updateInternals(){}; + UInt getUpdateRate(){return update_rate;}; protected: + NonLocalManager & manager; Real R; Real R2; + UInt update_rate; }; -class NonLocalNeighborhood : public NonLocalNeighborhoodBase { - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ -public: +class TestDamageWeightFunction : public TestWeightFunction { - NonLocalNeighborhood(NonLocalManager & manager, Real radius, const ID & weight_type, - const ElementTypeMapReal & quad_coordinates, - const ID & id = "neighborhood", - const MemoryID & memory_id = 0); - virtual ~NonLocalNeighborhood(); - - /* ------------------------------------------------------------------------ */ - /* Methods */ - /* ------------------------------------------------------------------------ */ public: + TestDamageWeightFunction(NonLocalManager & manager) : TestWeightFunction(manager){ + this->registerParam("damage_limit", this->damage_limit, 1., _pat_parsable, "Damage Threshold"); + this->setInternals();} - /// compute the weights for non-local averaging - void computeWeights(); - - /// save the pair of weights in a file - void saveWeights(const std::string & filename) const; +Real operator()(Real r, + const __attribute__((unused)) QuadraturePoint & q1, + const __attribute__((unused)) QuadraturePoint & q2) { + /// compute the weight + UInt quad = q2.global_num; - /// compute the non-local counter part for a given element type map - void weightedAvergageOnNeighbours(const ElementTypeMapReal & to_accumulate, - ElementTypeMapReal & accumulated, - UInt nb_degree_of_freedom, - const GhostType & ghost_type2) const; + if(q1 == q2) return 1.; - /* ------------------------------------------------------------------------ */ - /* Class Members */ - /* ------------------------------------------------------------------------ */ -private: - - /// Pointer to non-local manager class - NonLocalManager * non_local_manager; + Array & dam = (*(this->damage))(q2.type, q2.ghost_type); + Real D = dam(quad); + Real w = 0.; + if(D < damage_limit) { + Real alpha = std::max(0., 1. - r*r / this->R2); + w = alpha * alpha; + } + return w; +} + void updateInternals(){}; + UInt getUpdateRate(){return update_rate;}; - /// name of the type of weight function - const ID type; + void setInternals() { + this->damage = &(this->manager.registerWeightFunctionInternal("damage"));} - /// the weights associated to the pairs - Array * pair_weight[2]; +protected: - /// count the number of calls of computeStress - UInt compute_stress_calls; + /// limit at which a point is considered as complitely broken + Real damage_limit; - /// weight function - WeightFunction * weight_function; + /// internal pointer to the current damage vector + ElementTypeMapReal * damage; }; __END_AKANTU__ +/* -------------------------------------------------------------------------- */ +/* Implementation of template functions */ +/* -------------------------------------------------------------------------- */ +#include "non_local_neighborhood_tmpl.hh" +/* -------------------------------------------------------------------------- */ +/* inline functions */ +/* -------------------------------------------------------------------------- */ +#include "non_local_neighborhood_inline_impl.cc" + + + #endif /* __AKANTU_NON_LOCAL_NEIGHBORHOOD_HH__ */ + diff --git a/src/model/common/non_local_toolbox/non_local_neighborhood_base.cc b/src/model/common/non_local_toolbox/non_local_neighborhood_base.cc index 20387b491..a574fd7eb 100644 --- a/src/model/common/non_local_toolbox/non_local_neighborhood_base.cc +++ b/src/model/common/non_local_toolbox/non_local_neighborhood_base.cc @@ -1,278 +1,276 @@ /** * @file non_local_neighborhood_base.cc * @author Aurelia Isabel Cuba Ramos * @date Mon Sep 21 18:10:49 2015 * * @brief Implementation of non-local neighborhood base * * @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 "non_local_neighborhood_base.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ -NonLocalNeighborhoodBase::NonLocalNeighborhoodBase(const SolidMechanicsModel & model, - Real radius, +NonLocalNeighborhoodBase::NonLocalNeighborhoodBase(const SolidMechanicsModel & model, const ElementTypeMapReal & quad_coordinates, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), + Parsable(_st_non_local, id), model(model), - non_local_radius(radius), + non_local_radius(0.), spatial_grid(NULL), is_creating_grid(false), grid_synchronizer(NULL), - quad_coordinates(quad_coordinates){ + quad_coordinates(quad_coordinates), + spatial_dimension(this->model.getSpatialDimension()){ AKANTU_DEBUG_IN(); + this->registerParam("radius" , non_local_radius , 100., + _pat_parsable | _pat_readable , "Non local radius"); - this->initNeighborhood(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ NonLocalNeighborhoodBase::~NonLocalNeighborhoodBase() { AKANTU_DEBUG_IN(); delete spatial_grid; delete grid_synchronizer; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NonLocalNeighborhoodBase::initNeighborhood() { AKANTU_DEBUG_IN(); // Material::initMaterial(); Mesh & mesh = this->model.getMesh(); ElementTypeMap nb_ghost_protected; - UInt spatial_dimension = this->model.getSpatialDimension(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _ghost); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _ghost); for(; it != last_type; ++it) nb_ghost_protected(mesh.getNbElement(*it, _ghost), *it, _ghost); AKANTU_DEBUG_INFO("Creating the grid"); this->createGrid(); /// todo: create pairs, set radius of weight function, intialize weight function, compute weights AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------- */ void NonLocalNeighborhoodBase::createGrid() { AKANTU_DEBUG_IN(); const Real safety_factor = 1.2; // for the cell grid spacing Mesh & mesh = this->model.getMesh(); mesh.computeBoundingBox(); const Vector & lower_bounds = mesh.getLocalLowerBounds(); const Vector & upper_bounds = mesh.getLocalUpperBounds(); Vector center = 0.5 * (upper_bounds + lower_bounds); - UInt spatial_dimension = this->model.getSpatialDimension(); Vector spacing(spatial_dimension, this->non_local_radius * safety_factor); spatial_grid = new SpatialGrid(spatial_dimension, spacing, center); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NonLocalNeighborhoodBase::updatePairList() { AKANTU_DEBUG_IN(); - UInt spatial_dimension = this->model.getSpatialDimension(); //// loop over all quads -> all cells SpatialGrid::cells_iterator cell_it = spatial_grid->beginCells(); SpatialGrid::cells_iterator cell_end = spatial_grid->endCells(); Vector q1_coords(spatial_dimension); Vector q2_coords(spatial_dimension); QuadraturePoint q1; QuadraturePoint q2; UInt counter = 0; for (; cell_it != cell_end; ++cell_it) { AKANTU_DEBUG_INFO("Looping on next cell"); SpatialGrid::Cell::iterator first_quad = spatial_grid->beginCell(*cell_it); SpatialGrid::Cell::iterator last_quad = spatial_grid->endCell(*cell_it); for (;first_quad != last_quad; ++first_quad, ++counter){ q1 = *first_quad; Array::const_vector_iterator coords_type_1_it = this->quad_coordinates(q1.type, q1.ghost_type).begin(spatial_dimension); q1_coords = coords_type_1_it[q1.global_num]; AKANTU_DEBUG_INFO("Current quadrature point in this cell: " << q1); SpatialGrid::CellID cell_id = spatial_grid->getCellID(q1_coords); /// loop over all the neighbouring cells of the current quad SpatialGrid::neighbor_cells_iterator first_neigh_cell = spatial_grid->beginNeighborCells(cell_id); SpatialGrid::neighbor_cells_iterator last_neigh_cell = spatial_grid->endNeighborCells(cell_id); for (; first_neigh_cell != last_neigh_cell; ++first_neigh_cell) { SpatialGrid::Cell::iterator first_neigh_quad = spatial_grid->beginCell(*first_neigh_cell); SpatialGrid::Cell::iterator last_neigh_quad = spatial_grid->endCell(*first_neigh_cell); // loop over the quadrature point in the current neighboring cell for (;first_neigh_quad != last_neigh_quad; ++first_neigh_quad){ q2 = *first_neigh_quad; Array::const_vector_iterator coords_type_2_it = this->quad_coordinates(q2.type, q2.ghost_type).begin(spatial_dimension); q2_coords = coords_type_2_it[q2.global_num]; Real distance = q1_coords.distance(q2_coords); if(distance <= this->non_local_radius + Math::getTolerance() && (q2.ghost_type == _ghost || (q2.ghost_type == _not_ghost && q1.global_num <= q2.global_num))) { // storing only half lists pair_list[q2.ghost_type].push_back(std::make_pair(q1, q2)); } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NonLocalNeighborhoodBase::createGridSynchronizer() { this->is_creating_grid = true; std::set tags; tags.insert(_gst_mnl_for_average); tags.insert(_gst_mnl_weight); tags.insert(_gst_material_id); SynchronizerRegistry & synch_registry = this->model.getSynchronizerRegistry(); std::stringstream sstr; sstr << getID() << ":grid_synchronizer"; this->grid_synchronizer = GridSynchronizer::createGridSynchronizer(this->model.getMesh(), *spatial_grid, sstr.str(), &synch_registry, tags); this->is_creating_grid = false; } /* -------------------------------------------------------------------------- */ void NonLocalNeighborhoodBase::savePairs(const std::string & filename) const { std::ofstream pout; std::stringstream sstr; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int prank = comm.whoAmI(); sstr << filename << "." << prank; pout.open(sstr.str().c_str()); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; PairList::const_iterator first_pair = pair_list[ghost_type2].begin(); PairList::const_iterator last_pair = pair_list[ghost_type2].end(); for(;first_pair != last_pair; ++first_pair) { const QuadraturePoint & q1 = first_pair->first; const QuadraturePoint & q2 = first_pair->second; pout << q1 << " " << q2 << " " << std::endl; } } } /* -------------------------------------------------------------------------- */ void NonLocalNeighborhoodBase::saveNeighborCoords(const std::string & filename) const { /// this function is not optimazed and only used for tests on small meshes /// @todo maybe optimize this function for better performance? - UInt spatial_dimension = this->model.getSpatialDimension(); Vector q1_coords(spatial_dimension); Vector q2_coords(spatial_dimension); QuadraturePoint q1; QuadraturePoint q2; std::ofstream pout; std::stringstream sstr; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int prank = comm.whoAmI(); sstr << filename << "." << prank; pout.open(sstr.str().c_str()); /// loop over all the quads and write the position of their neighbors SpatialGrid::cells_iterator cell_it = spatial_grid->beginCells(); SpatialGrid::cells_iterator cell_end = spatial_grid->endCells(); for (; cell_it != cell_end; ++cell_it) { SpatialGrid::Cell::iterator first_quad = spatial_grid->beginCell(*cell_it); SpatialGrid::Cell::iterator last_quad = spatial_grid->endCell(*cell_it); for (;first_quad != last_quad; ++first_quad){ q1 = *first_quad; Array::const_vector_iterator coords_type_1_it = this->quad_coordinates(q1.type, q1.ghost_type).begin(spatial_dimension); q1_coords = coords_type_1_it[q1.global_num]; pout << "#neighbors for quad " << q1.global_num << std::endl; pout << q1_coords << std::endl; for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; PairList::const_iterator first_pair = pair_list[ghost_type2].begin(); PairList::const_iterator last_pair = pair_list[ghost_type2].end(); for(;first_pair != last_pair; ++first_pair) { if (q1 == first_pair->first && first_pair->second != q1) { q2 = first_pair->second; Array::const_vector_iterator coords_type_2_it = this->quad_coordinates(q2.type, q2.ghost_type).begin(spatial_dimension); q2_coords = coords_type_2_it[q2.global_num]; pout << q2_coords << std::endl; } if (q1 == first_pair->second && first_pair->first != q1) { q2 = first_pair->first; Array::const_vector_iterator coords_type_2_it = this->quad_coordinates(q2.type, q2.ghost_type).begin(spatial_dimension); q2_coords = coords_type_2_it[q2.global_num]; pout << q2_coords << std::endl; } } } } } } __END_AKANTU__ diff --git a/src/model/common/non_local_toolbox/non_local_neighborhood_base.hh b/src/model/common/non_local_toolbox/non_local_neighborhood_base.hh index 364679be7..43dd78da3 100644 --- a/src/model/common/non_local_toolbox/non_local_neighborhood_base.hh +++ b/src/model/common/non_local_toolbox/non_local_neighborhood_base.hh @@ -1,149 +1,167 @@ /** * @file non_local_neighborhood_base.hh * @author Aurelia Isabel Cuba Ramos * @date Mon Sep 21 15:43:26 2015 * * @brief Non-local neighborhood base class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_NON_LOCAL_NEIGHBORHOOD_BASE_HH__ #define __AKANTU_NON_LOCAL_NEIGHBORHOOD_BASE_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "solid_mechanics_model.hh" #include "aka_grid_dynamic.hh" #include "grid_synchronizer.hh" #include "aka_memory.hh" +#include "data_accessor.hh" +#include "parsable.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ -class NonLocalNeighborhoodBase : public Memory { +class NonLocalNeighborhoodBase : public Memory, + public DataAccessor, + public Parsable{ /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - NonLocalNeighborhoodBase(const SolidMechanicsModel & model, Real radius, + NonLocalNeighborhoodBase(const SolidMechanicsModel & model, const ElementTypeMapReal & quad_coordinates, const ID & id = "neighborhood", const MemoryID & memory_id = 0); virtual ~NonLocalNeighborhoodBase(); typedef std::vector< std::pair > PairList; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: + + /// intialize the neighborhood + void initNeighborhood(); + /// initialize the material computed parameter inline void insertQuad(const QuadraturePoint & quad, const Vector & coords); /// create the pairs of quadrature points void updatePairList(); /// save the pairs of quadrature points in a file void savePairs(const std::string & filename) const; /// save the coordinates of all neighbors of a quad void saveNeighborCoords(const std::string & filename) const; /// create grid synchronizer and exchange ghost cells void createGridSynchronizer(); /// compute weights, for instance needed for non-local damage computation virtual void computeWeights() {}; /// compute the non-local counter part for a given element type map - void weightedAvergageOnNeighbours(const ElementTypeMapReal & to_accumulate, - ElementTypeMapReal & accumulated, - UInt nb_degree_of_freedom, - const GhostType & ghost_type2) const {}; + virtual void weightedAverageOnNeighbours(const ElementTypeMapReal & to_accumulate, + ElementTypeMapReal & accumulated, + UInt nb_degree_of_freedom, + const GhostType & ghost_type2) const {}; -protected: + /// update the weights for the non-local averaging + virtual void updateWeights() {}; - /// intialize the neighborhood - void initNeighborhood(); +protected: /// create the grid void createGrid(); void cleanupExtraGhostElement(const ElementTypeMap & nb_ghost_protected) {}; - // virtual inline UInt getNbDataForElements(const Array & elements, - // SynchronizationTag tag) const; + virtual inline UInt getNbDataForElements(const Array & elements, + SynchronizationTag tag) const {return 0; }; - // virtual inline void packElementData(CommunicationBuffer & buffer, - // const Array & elements, - // SynchronizationTag tag) const; + virtual inline void packElementData(CommunicationBuffer & buffer, + const Array & elements, + SynchronizationTag tag) const {}; - // virtual inline void unpackElementData(CommunicationBuffer & buffer, - // const Array & elements, - // SynchronizationTag tag); + virtual inline void unpackElementData(CommunicationBuffer & buffer, + const Array & elements, + SynchronizationTag tag) {}; // virtual inline void onElementsAdded(const Array & element_list); // virtual inline void onElementsRemoved(const Array & element_list, // const ElementTypeMapArray & new_numbering, // const RemovedElementsEvent & event) {}; +/* -------------------------------------------------------------------------- */ +/* Accessors */ +/* -------------------------------------------------------------------------- */ +public: + AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); + AKANTU_GET_MACRO(Model, model, const SolidMechanicsModel &); + /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// the model to which the neighborhood belongs const SolidMechanicsModel & model; /// Radius of non-local neighborhood Real non_local_radius; /** * the pairs of quadrature points * 0: not ghost to not ghost * 1: not ghost to ghost */ PairList pair_list[2]; /// the regular grid to construct/update the pair lists SpatialGrid * spatial_grid; bool is_creating_grid; /// the grid synchronizer for parallel computations GridSynchronizer * grid_synchronizer; /// the quadrature point positions const ElementTypeMapReal & quad_coordinates; + + /// the spatial dimension of the problem + const UInt spatial_dimension; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "non_local_neighborhood_base_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_NON_LOCAL_NEIGHBORHOOD_BASE_HH__ */ diff --git a/src/model/common/non_local_toolbox/non_local_neighborhood_base_inline_impl.cc b/src/model/common/non_local_toolbox/non_local_neighborhood_base_inline_impl.cc index 95e5f84fd..a2bf76460 100644 --- a/src/model/common/non_local_toolbox/non_local_neighborhood_base_inline_impl.cc +++ b/src/model/common/non_local_toolbox/non_local_neighborhood_base_inline_impl.cc @@ -1,46 +1,45 @@ /** * @file non_local_neighborhood_base_inline_impl.cc * @author Aurelia Isabel Cuba Ramos * @date Mon Sep 21 18:25:31 2015 * * @brief Inline implementation of Non-local neighborhood base functions * * @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__ //// includes /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline void NonLocalNeighborhoodBase::insertQuad(const QuadraturePoint & quad, const Vector & coords) { - Vector test = quad.getPosition(); this->spatial_grid->insert(quad, coords); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ diff --git a/src/model/common/non_local_toolbox/non_local_neighborhood_inline_impl.cc b/src/model/common/non_local_toolbox/non_local_neighborhood_inline_impl.cc new file mode 100644 index 000000000..9c81afc57 --- /dev/null +++ b/src/model/common/non_local_toolbox/non_local_neighborhood_inline_impl.cc @@ -0,0 +1,95 @@ +/** + * @file non_local_neighborhood_inline_impl.cc + * @author Aurelia Isabel Cuba Ramos + * @date Mon Oct 5 09:36:33 2015 + * + * @brief Implementation of inline functions of non-local neighborhood class + * + * @section LICENSE + * + * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + */ + +/* -------------------------------------------------------------------------- */ +#ifndef __AKANTU_NON_LOCAL_NEIGHBORHOOD_INLINE_IMPL_CC__ +#define __AKANTU_NON_LOCAL_NEIGHBORHOOD_INLINE_IMPL_CC__ + +__BEGIN_AKANTU__ +/* -------------------------------------------------------------------------- */ +template +inline UInt NonLocalNeighborhood::getNbDataForElements(const Array & elements, + SynchronizationTag tag) const { + // UInt nb_quadrature_points = this->getModel().getNbQuadraturePoints(elements); + // UInt size = 0; + + // if(tag == _gst_mnl_for_average) { + // typename std::map::const_iterator it = non_local_variables.begin(); + // typename std::map::const_iterator end = non_local_variables.end(); + + // for(;it != end; ++it) { + // const NonLocalVariable & non_local_variable = it->second; + // size += non_local_variable.nb_component * sizeof(Real) * nb_quadrature_points; + // } + // } + + // size += weight_func->getNbDataForElements(elements, tag); + + return 0; +} + +/* -------------------------------------------------------------------------- */ +template +inline void NonLocalNeighborhood::packElementData(CommunicationBuffer & buffer, + const Array & elements, + SynchronizationTag tag) const { + // if(tag == _gst_mnl_for_average) { + // typename std::map::const_iterator it = non_local_variables.begin(); + // typename std::map::const_iterator end = non_local_variables.end(); + + // for(;it != end; ++it) { + // const NonLocalVariable & non_local_variable = it->second; + // this->packElementDataHelper(*non_local_variable.local, + // buffer, elements); + // } + // } + + // weight_func->packElementData(buffer, elements, tag); +} + +/* -------------------------------------------------------------------------- */ +template +inline void NonLocalNeighborhood::unpackElementData(CommunicationBuffer & buffer, + const Array & elements, + SynchronizationTag tag) { + // if(tag == _gst_mnl_for_average) { + // typename std::map::iterator it = non_local_variables.begin(); + // typename std::map::iterator end = non_local_variables.end(); + + // for(;it != end; ++it) { + // NonLocalVariable & non_local_variable = it->second; + // this->unpackElementDataHelper(*non_local_variable.local, + // buffer, elements); + // } + // } + + // weight_func->unpackElementData(buffer, elements, tag); +} + +__END_AKANTU__ + +#endif /* __AKANTU_NON_LOCAL_NEIGHBORHOOD_INLINE_IMPL_CC__ */ diff --git a/src/model/common/non_local_toolbox/non_local_neighborhood.cc b/src/model/common/non_local_toolbox/non_local_neighborhood_tmpl.hh similarity index 68% rename from src/model/common/non_local_toolbox/non_local_neighborhood.cc rename to src/model/common/non_local_toolbox/non_local_neighborhood_tmpl.hh index 40f5e9808..7f8f9d2d8 100644 --- a/src/model/common/non_local_toolbox/non_local_neighborhood.cc +++ b/src/model/common/non_local_toolbox/non_local_neighborhood_tmpl.hh @@ -1,260 +1,272 @@ /** - * @file non_local_neighborhood.cc + * @file non_local_neighborhood_tmpl.hh * @author Aurelia Isabel Cuba Ramos * @date Sat Sep 26 18:43:30 2015 * * @brief Implementation of class non-local neighborhood * * @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 "non_local_neighborhood.hh" #include "non_local_manager.hh" /* -------------------------------------------------------------------------- */ +#ifndef __AKANTU_NON_LOCAL_NEIGHBORHOOD_TMPL_HH__ +#define __AKANTU_NON_LOCAL_NEIGHBORHOOD_TMPL_HH__ + __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ -NonLocalNeighborhood::NonLocalNeighborhood(NonLocalManager & manager, - Real radius, - const ID & type, - const ElementTypeMapReal & quad_coordinates, - const ID & id, - const MemoryID & memory_id) : - NonLocalNeighborhoodBase(manager.getModel(), radius, quad_coordinates, id, memory_id), - non_local_manager(&manager), - type(type) { +template +NonLocalNeighborhood::NonLocalNeighborhood(NonLocalManager & manager, + const ElementTypeMapReal & quad_coordinates, + const ID & id, + const MemoryID & memory_id) : + NonLocalNeighborhoodBase(manager.getModel(), quad_coordinates, id, memory_id), + non_local_manager(&manager) { AKANTU_DEBUG_IN(); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type = (GhostType) gt; pair_weight[ghost_type] = NULL; } - // if (weight_func_type == "base_wf") - this->weight_function = new WeightFunction(); - this->weight_function->setRadius(radius); - - // case damage_wf: - // this->weight_function = new DamagedWeightFunction(); break; - // case remove_wf: - // this->weight_function = new RemoveDamagedWeightFunction(); break; - // case stress_wf: - // this->weight_function = new StressBasedWeightFunction(); break + this->weight_function = new WeightFunction(this->getNonLocalManager()); + this->registerSubSection(_st_weight_function, "weight_parameter", *weight_function); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -NonLocalNeighborhood::~NonLocalNeighborhood() { +template +NonLocalNeighborhood::~NonLocalNeighborhood() { AKANTU_DEBUG_IN(); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type = (GhostType) gt; delete pair_weight[ghost_type]; } delete weight_function; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void NonLocalNeighborhood::computeWeights() { +template +void NonLocalNeighborhood::computeWeights() { AKANTU_DEBUG_IN(); - UInt spatial_dimension = this->model.getSpatialDimension(); - Vector q1_coord(spatial_dimension); - Vector q2_coord(spatial_dimension); + this->weight_function->setRadius(this->non_local_radius); + Vector q1_coord(this->spatial_dimension); + Vector q2_coord(this->spatial_dimension); UInt nb_weights_per_pair = 2; /// w1: q1->q2, w2: q2->q1 /// get the elementtypemap for the neighborhood volume for each quadrature point ElementTypeMapReal & quadrature_points_volumes = this->non_local_manager->getVolumes(); /// update the internals of the weight function if applicable (not /// all the weight functions have internals and do noting in that /// case) weight_function->updateInternals(); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type = (GhostType) gt; /// allocate the array to store the weight, if it doesn't exist already if(!(pair_weight[ghost_type])) { std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; std::stringstream sstr; sstr << this->id <<":pair_weight:" << ghost_id; pair_weight[ghost_type] = new Array(0, nb_weights_per_pair, sstr.str()); } /// resize the array to the correct size pair_weight[ghost_type]->resize(pair_list[ghost_type].size()); /// set entries to zero pair_weight[ghost_type]->clear(); /// loop over all pairs in the current pair list array and their corresponding weights PairList::const_iterator first_pair = pair_list[ghost_type].begin(); PairList::const_iterator last_pair = pair_list[ghost_type].end(); Array::vector_iterator weight_it = pair_weight[ghost_type]->begin(nb_weights_per_pair); // Compute the weights for(;first_pair != last_pair; ++first_pair, ++weight_it) { Vector & weight = *weight_it; const QuadraturePoint & q1 = first_pair->first; const QuadraturePoint & q2 = first_pair->second; /// get the coordinates for the given pair of quads - Array::const_vector_iterator coords_type_1_it = this->quad_coordinates(q1.type, q1.ghost_type).begin(spatial_dimension); + Array::const_vector_iterator coords_type_1_it = this->quad_coordinates(q1.type, q1.ghost_type).begin(this->spatial_dimension); q1_coord = coords_type_1_it[q1.global_num]; - Array::const_vector_iterator coords_type_2_it = this->quad_coordinates(q2.type, q2.ghost_type).begin(spatial_dimension); + Array::const_vector_iterator coords_type_2_it = this->quad_coordinates(q2.type, q2.ghost_type).begin(this->spatial_dimension); q2_coord = coords_type_2_it[q2.global_num]; Array & quad_volumes_1 = quadrature_points_volumes(q1.type, q1.ghost_type); const Array & jacobians_2 = - this->non_local_manager->getJacobians(q2.type, q2.ghost_type); + this->non_local_manager->getJacobians(q2.type, q2.ghost_type); const Real & q2_wJ = jacobians_2(q2.global_num); /// compute distance between the two quadrature points Real r = q1_coord.distance(q2_coord); /// compute the weight for averaging on q1 based on the distance Real w1 = this->weight_function->operator()(r, q1, q2); weight(0) = q2_wJ * w1; quad_volumes_1(q1.global_num) += weight(0); if(q2.ghost_type != _ghost && q1.global_num != q2.global_num) { - const Array & jacobians_1 = - this->non_local_manager->getJacobians(q1.type, q1.ghost_type); - Array & quad_volumes_2 = - quadrature_points_volumes(q2.type, q2.ghost_type); + const Array & jacobians_1 = + this->non_local_manager->getJacobians(q1.type, q1.ghost_type); + Array & quad_volumes_2 = + quadrature_points_volumes(q2.type, q2.ghost_type); /// compute the weight for averaging on q2 - const Real & q1_wJ = jacobians_1(q1.global_num); - Real w2 = this->weight_function->operator()(r, q2, q1); - weight(1) = q1_wJ * w2; - quad_volumes_2(q2.global_num) += weight(1); + const Real & q1_wJ = jacobians_1(q1.global_num); + Real w2 = this->weight_function->operator()(r, q2, q1); + weight(1) = q1_wJ * w2; + quad_volumes_2(q2.global_num) += weight(1); } else - weight(1) = 0.; + weight(1) = 0.; } } /// normalize the weights for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; PairList::const_iterator first_pair = pair_list[ghost_type2].begin(); PairList::const_iterator last_pair = pair_list[ghost_type2].end(); Array::vector_iterator weight_it = pair_weight[ghost_type2]->begin(nb_weights_per_pair); // Compute the weights for(;first_pair != last_pair; ++first_pair, ++weight_it) { Vector & weight = *weight_it; const QuadraturePoint & q1 = first_pair->first; const QuadraturePoint & q2 = first_pair->second; Array & quad_volumes_1 = quadrature_points_volumes(q1.type, q1.ghost_type); Array & quad_volumes_2 = quadrature_points_volumes(q2.type, q2.ghost_type); Real q1_volume = quad_volumes_1(q1.global_num); weight(0) *= 1. / q1_volume; if(ghost_type2 != _ghost) { - Real q2_volume = quad_volumes_2(q2.global_num); - weight(1) *= 1. / q2_volume; + Real q2_volume = quad_volumes_2(q2.global_num); + weight(1) *= 1. / q2_volume; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void NonLocalNeighborhood::saveWeights(const std::string & filename) const { +template +void NonLocalNeighborhood::saveWeights(const std::string & filename) const { std::ofstream pout; std::stringstream sstr; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int prank = comm.whoAmI(); sstr << filename << "." << prank; pout.open(sstr.str().c_str()); for (UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type = (GhostType) gt; - AKANTU_DEBUG_ASSERT((pair_weight[ghost_type]), "the weights have not been computed yet"); + AKANTU_DEBUG_ASSERT((pair_weight[ghost_type]), "the weights have not been computed yet"); - Array & weights = *(pair_weight[ghost_type]); - Array::const_vector_iterator weights_it = weights.begin(2); - for (UInt i = 0; i < weights.getSize(); ++i, ++weights_it) - pout << "w1: " << (*weights_it)(0) <<" w2: " << (*weights_it)(1) << std::endl; + Array & weights = *(pair_weight[ghost_type]); + Array::const_vector_iterator weights_it = weights.begin(2); + for (UInt i = 0; i < weights.getSize(); ++i, ++weights_it) + pout << "w1: " << (*weights_it)(0) <<" w2: " << (*weights_it)(1) << std::endl; } } /* -------------------------------------------------------------------------- */ -void NonLocalNeighborhood::weightedAvergageOnNeighbours(const ElementTypeMapReal & to_accumulate, - ElementTypeMapReal & accumulated, - UInt nb_degree_of_freedom, - const GhostType & ghost_type2) const { +template +void NonLocalNeighborhood::weightedAverageOnNeighbours(const ElementTypeMapReal & to_accumulate, + ElementTypeMapReal & accumulated, + UInt nb_degree_of_freedom, + const GhostType & ghost_type2) const { AKANTU_DEBUG_IN(); PairList::const_iterator first_pair = pair_list[ghost_type2].begin(); PairList::const_iterator last_pair = pair_list[ghost_type2].end(); Array::vector_iterator weight_it = pair_weight[ghost_type2]->begin(2); // Compute the weights for(;first_pair != last_pair; ++first_pair, ++weight_it) { Vector & weight = *weight_it; const QuadraturePoint & q1 = first_pair->first; const QuadraturePoint & q2 = first_pair->second; const Array & to_acc_1 = to_accumulate(q1.type, q1.ghost_type); Array & acc_1 = accumulated(q1.type, q1.ghost_type); const Array & to_acc_2 = to_accumulate(q2.type, q2.ghost_type); Array & acc_2 = accumulated(q2.type, q2.ghost_type); for(UInt d = 0; d < nb_degree_of_freedom; ++d) { acc_1(q1.global_num, d) += weight(0) * to_acc_2(q2.global_num, d); } if(ghost_type2 != _ghost) { for(UInt d = 0; d < nb_degree_of_freedom; ++d) { acc_2(q2.global_num, d) += weight(1) * to_acc_1(q1.global_num, d); } } } AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +template +void NonLocalNeighborhood::updateWeights() { + // Update the weights for the non local variable averaging + if(this->weight_function->getUpdateRate() && + (this->non_local_manager->getNbStressCalls() % this->weight_function->getUpdateRate() == 0)) { + // this->model->getSynchronizerRegistry().asynchronousSynchronize(_gst_mnl_weight); + // this->model->getSynchronizerRegistry().waitEndSynchronize(_gst_mnl_weight); + this->computeWeights(); + } +} + +/* -------------------------------------------------------------------------- */ + __END_AKANTU__ + +#endif /* __AKANTU_NON_LOCAL_NEIGHBORHOOD_TMPL__ */ diff --git a/src/model/solid_mechanics/materials/material_non_local.hh b/src/model/solid_mechanics/materials/material_non_local.hh index 0396735c9..1d59de45c 100644 --- a/src/model/solid_mechanics/materials/material_non_local.hh +++ b/src/model/solid_mechanics/materials/material_non_local.hh @@ -1,189 +1,194 @@ /** * @file material_non_local.hh * * @author Nicolas Richart * * @date creation: Wed Aug 31 2011 * @date last modification: Thu Jun 05 2014 * * @brief Material class that handle the non locality of a law for example damage. * * @section LICENSE * * Copyright (©) 2014 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 "material.hh" #include "aka_grid_dynamic.hh" #include "fe_engine.hh" #include "base_weight_function.hh" namespace akantu { class GridSynchronizer; } /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_NON_LOCAL_HH__ #define __AKANTU_MATERIAL_NON_LOCAL_HH__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template class MaterialNonLocal : public virtual Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialNonLocal(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialNonLocal(); typedef std::vector< std::pair > PairList; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the material computed parameter virtual void initMaterial(); virtual void updateResidual(GhostType ghost_type); virtual void computeAllNonLocalStresses(GhostType ghost_type = _not_ghost); void savePairs(const std::string & filename) const; void saveWeights(const std::string & filename) const; void neighbourhoodStatistics(const std::string & filename) const; + /// insert the quadrature points in the neighborhoods of the non-local manager virtual void insertQuadsInNeighborhoods(GhostType ghost_type = _not_ghost); + /// update the values in the non-local internal fields + void updateNonLocalInternals(ElementTypeMapReal & non_local_flattened, const ID & field_id, + const UInt nb_component); + + template + void weightedAvergageOnNeighbours(const InternalField & to_accumulate, + InternalField & accumulated, + UInt nb_degree_of_freedom, + GhostType ghost_type2 = _not_ghost) const; + + /// constitutive law + virtual void computeNonLocalStresses(GhostType ghost_type = _not_ghost) = 0; + protected: void updatePairList(const ElementTypeMapArray & quadrature_points_coordinates); void computeWeights(const ElementTypeMapArray & quadrature_points_coordinates); void createCellList(ElementTypeMapArray & quadrature_points_coordinates); void cleanupExtraGhostElement(const ElementTypeMap & nb_ghost_protected); void fillCellList(const ElementTypeMapArray & quadrature_points_coordinates, const GhostType & ghost_type); - /// constitutive law - virtual void computeNonLocalStresses(GhostType ghost_type = _not_ghost) = 0; - - template - void weightedAvergageOnNeighbours(const InternalField & to_accumulate, - InternalField & accumulated, - UInt nb_degree_of_freedom, - GhostType ghost_type2 = _not_ghost) const; - virtual inline UInt getNbDataForElements(const Array & elements, SynchronizationTag tag) const; virtual inline void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const; virtual inline void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag); virtual inline void onElementsAdded(const Array & element_list); virtual inline void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: void registerNonLocalVariable(InternalField & local, InternalField & non_local, UInt nb_degree_of_freedom) { ID id = local.getID(); NonLocalVariable & non_local_variable = non_local_variables[id]; non_local_variable.local = &local; non_local_variable.non_local = &non_local; non_local_variable.nb_component = nb_degree_of_freedom; } AKANTU_GET_MACRO(PairList, pair_list, const PairList &) Real getRadius() const { return weight_func->getRadius(); } AKANTU_GET_MACRO(CellList, *spatial_grid, const SpatialGrid &) /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// the weight function used WeightFunction * weight_func; private: /** * the pairs of quadrature points * 0: not ghost to not ghost * 1: not ghost to ghost */ PairList pair_list[2]; /// the weights associated to the pairs Array * pair_weight[2]; /// the regular grid to construct/update the pair lists SpatialGrid * spatial_grid; /// the types of the existing pairs typedef std::set< std::pair > pair_type; pair_type existing_pairs[2]; /// count the number of calls of computeStress UInt compute_stress_calls; struct NonLocalVariable { InternalField * local; InternalField * non_local; UInt nb_component; }; std::map non_local_variables; bool is_creating_grid; GridSynchronizer * grid_synchronizer; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_non_local_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_NON_LOCAL_HH__ */ 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 580297832..bd4dbe500 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,908 +1,937 @@ /** * @file material_non_local_inline_impl.cc * * @author Nicolas Richart * * @date creation: Wed Aug 31 2011 * @date last modification: Mon Jun 23 2014 * * @brief Non-local inline implementation * * @section LICENSE * * Copyright (©) 2014 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 "integrator.hh" #include "dumper_paraview.hh" #include "non_local_manager.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #if defined(AKANTU_DEBUG_TOOLS) # include "aka_debug_tools.hh" #endif __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialNonLocal::MaterialNonLocal(SolidMechanicsModel & model, const ID & id) : Material(model, id), weight_func(NULL), spatial_grid(NULL), compute_stress_calls(0), is_creating_grid(false), grid_synchronizer(NULL) { AKANTU_DEBUG_IN(); this->model->getNonLocalManager().registerNonLocalMaterial(*this); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type = (GhostType) gt; pair_weight[ghost_type] = NULL; } this->is_non_local = true; this->weight_func = new WeightFunction(*this); this->registerSubSection(_st_non_local, "weight_function", *weight_func); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template MaterialNonLocal::~MaterialNonLocal() { AKANTU_DEBUG_IN(); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type = (GhostType) gt; delete pair_weight[ghost_type]; } delete spatial_grid; delete weight_func; delete grid_synchronizer; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::initMaterial() { AKANTU_DEBUG_IN(); this->insertQuadsInNeighborhoods(_not_ghost); // Material::initMaterial(); Mesh & mesh = this->model->getFEEngine().getMesh(); InternalField quadrature_points_coordinates("quadrature_points_coordinates_tmp_nl", *this); quadrature_points_coordinates.initialize(spatial_dimension); ElementTypeMap nb_ghost_protected; Mesh::type_iterator it = mesh.firstType(spatial_dimension, _ghost); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _ghost); for(; it != last_type; ++it) nb_ghost_protected(mesh.getNbElement(*it, _ghost), *it, _ghost); AKANTU_DEBUG_INFO("Creating cell list"); createCellList(quadrature_points_coordinates); AKANTU_DEBUG_INFO("Creating pairs"); updatePairList(quadrature_points_coordinates); #if not defined(AKANTU_NDEBUG) if(AKANTU_DEBUG_TEST(dblDump)) neighbourhoodStatistics("material_non_local.stats"); #endif AKANTU_DEBUG_INFO("Cleaning extra ghosts"); /// cleanupExtraGhostElement(nb_ghost_protected); AKANTU_DEBUG_INFO("Computing weights"); weight_func->setRadius(weight_func->getRadius()); weight_func->init(); computeWeights(quadrature_points_coordinates); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::cleanupExtraGhostElement(const ElementTypeMap & nb_ghost_protected) { AKANTU_DEBUG_IN(); // Create list of element to keep std::set relevant_ghost_element; PairList::const_iterator first_pair = pair_list[_ghost].begin(); PairList::const_iterator last_pair = pair_list[_ghost].end(); for(;first_pair != last_pair; ++first_pair) { const QuadraturePoint & q2 = first_pair->second; relevant_ghost_element.insert(q2); } // Create list of element to remove and new numbering for element to keep Mesh & mesh = this->model->getFEEngine().getMesh(); std::set ghost_to_erase; Mesh::type_iterator it = mesh.firstType(spatial_dimension, _ghost); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _ghost); RemovedElementsEvent remove_elem(mesh); Element element_local; // member element corresponds to global element number element_local.ghost_type = _ghost; for(; it != last_type; ++it) { element_local.type = *it; UInt nb_ghost_elem_global = mesh.getNbElement(*it, _ghost); UInt nb_ghost_elem_protected = 0; try { nb_ghost_elem_protected = nb_ghost_protected(*it, _ghost); } catch (...) {} if(!remove_elem.getNewNumbering().exists(*it, _ghost)) remove_elem.getNewNumbering().alloc(nb_ghost_elem_global, 1, *it, _ghost); else remove_elem.getNewNumbering(*it, _ghost).resize(nb_ghost_elem_global); Array & elem_filter = element_filter(*it, _ghost); UInt nb_ghost_elem_local = elem_filter.getSize(); Array & new_numbering = remove_elem.getNewNumbering(*it, _ghost); for (UInt g = 0; g < nb_ghost_elem_local; ++g) { element_local.element = g; Element element_global = this->convertToGlobalElement(element_local); if (element_global.element >= nb_ghost_elem_protected && relevant_ghost_element.find(element_local) == relevant_ghost_element.end()) { // (std::find(relevant_ghost_element.begin(), // relevant_ghost_element.end(), // element_local) == relevant_ghost_element.end())) { std::cout<< "element removed" << std::endl; remove_elem.getList().push_back(element_global); new_numbering(element_global.element) = UInt(-1); } } UInt ng = 0; for (UInt g = 0; g < nb_ghost_elem_global; ++g) { if (new_numbering(g) != UInt(-1)) { new_numbering(g) = ng; ++ng; } } } mesh.sendEvent(remove_elem); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::createCellList(ElementTypeMapArray & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); const Real safety_factor = 1.2; // for the cell grid spacing Mesh & mesh = this->model->getFEEngine().getMesh(); mesh.computeBoundingBox(); const Vector & lower_bounds = mesh.getLocalLowerBounds(); const Vector & upper_bounds = mesh.getLocalUpperBounds(); Vector center = 0.5 * (upper_bounds + lower_bounds); Vector spacing(spatial_dimension, weight_func->getRadius() * safety_factor); spatial_grid = new SpatialGrid(spatial_dimension, spacing, center); this->fem->computeQuadraturePointsCoordinates(quadrature_points_coordinates, &element_filter); this->fillCellList(quadrature_points_coordinates, _not_ghost); is_creating_grid = true; std::set tags; tags.insert(_gst_mnl_for_average); tags.insert(_gst_mnl_weight); tags.insert(_gst_material_id); SynchronizerRegistry & synch_registry = this->model->getSynchronizerRegistry(); std::stringstream sstr; sstr << getID() << ":grid_synchronizer"; grid_synchronizer = GridSynchronizer::createGridSynchronizer(mesh, *spatial_grid, sstr.str(), &synch_registry, tags); is_creating_grid = false; this->fem->computeQuadraturePointsCoordinates(quadrature_points_coordinates, &element_filter); fillCellList(quadrature_points_coordinates, _ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::fillCellList(const ElementTypeMapArray & quadrature_points_coordinates, const GhostType & ghost_type) { QuadraturePoint q; q.ghost_type = ghost_type; q.kind = _ek_regular; Mesh::type_iterator it = this->element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = this->element_filter.lastType (spatial_dimension, ghost_type); for(; it != last_type; ++it) { q.type = *it; Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_quad = this->model->getFEEngine().getNbQuadraturePoints(*it, ghost_type); const Array & quads = quadrature_points_coordinates(*it, ghost_type); q.type = *it; Array::const_vector_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.global_num = q.element * nb_quad + nq; //q.setPosition(*quad); spatial_grid->insert(q, *quad); ++quad; } ++elem; } } } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::updatePairList(const ElementTypeMapArray & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); GhostType ghost_type = _not_ghost; QuadraturePoint q1; q1.ghost_type = ghost_type; // generate the pair of neighbor depending of the cell_list Mesh::type_iterator it = this->element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = this->element_filter.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { // Preparing datas const Array & quads = quadrature_points_coordinates(*it, ghost_type); Array::const_vector_iterator q1_coord_it = quads.begin(spatial_dimension); Array::const_vector_iterator last_quad = quads.end(spatial_dimension); q1.type = *it; q1.global_num = 0; // loop over quad points for(;q1_coord_it != last_quad; ++q1_coord_it, ++(q1.global_num)) { UInt nb_quad1 = this->model->getFEEngine().getNbQuadraturePoints(q1.type, q1.ghost_type); q1.element = q1.global_num / nb_quad1; q1.num_point = q1.global_num % nb_quad1; const Vector & q1_coord = *q1_coord_it; SpatialGrid::CellID cell_id = spatial_grid->getCellID(q1_coord); SpatialGrid::neighbor_cells_iterator first_neigh_cell = spatial_grid->beginNeighborCells(cell_id); SpatialGrid::neighbor_cells_iterator last_neigh_cell = spatial_grid->endNeighborCells(cell_id); // loop over neighbors cells of the one containing the current quadrature // point for (; first_neigh_cell != last_neigh_cell; ++first_neigh_cell) { SpatialGrid::Cell::iterator first_neigh_quad = spatial_grid->beginCell(*first_neigh_cell); SpatialGrid::Cell::iterator last_neigh_quad = spatial_grid->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 q2 = this->convertToLocalPoint(*first_neigh_quad); Array::const_vector_iterator q2_coord_it = quadrature_points_coordinates(q2.type, q2.ghost_type).begin(spatial_dimension); const Vector & q2_coord = q2_coord_it[q2.global_num]; Real distance = q1_coord.distance(q2_coord); if(distance <= weight_func->getRadius() && (q2.ghost_type == _ghost || (q2.ghost_type == _not_ghost && q1.global_num <= q2.global_num))) { // storing only half lists pair_list[q2.ghost_type].push_back(std::make_pair(q1, q2)); } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::computeWeights(const ElementTypeMapArray & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); InternalField quadrature_points_volumes("quadrature_points_volumes", *this); quadrature_points_volumes.initialize(1); const FEEngine & fem = this->model->getFEEngine(); weight_func->updateInternals(); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; if(!(pair_weight[ghost_type2])) { std::string ghost_id = ""; if (ghost_type2 == _ghost) ghost_id = ":ghost"; std::stringstream sstr; sstr << getID() << ":pair_weight:" << ghost_id; pair_weight[ghost_type2] = new Array(0, 2, sstr.str()); } pair_weight[ghost_type2]->resize(pair_list[ghost_type2].size()); pair_weight[ghost_type2]->clear(); PairList::const_iterator first_pair = pair_list[ghost_type2].begin(); PairList::const_iterator last_pair = pair_list[ghost_type2].end(); Array::vector_iterator weight_it = pair_weight[ghost_type2]->begin(2); // Compute the weights for(;first_pair != last_pair; ++first_pair, ++weight_it) { Vector & weight = *weight_it; const QuadraturePoint & lq1 = first_pair->first; const QuadraturePoint & lq2 = first_pair->second; QuadraturePoint gq1 = this->convertToGlobalPoint(lq1); QuadraturePoint gq2 = this->convertToGlobalPoint(lq2); // const Real q2_wJ = fem.getIntegratorInterface().getJacobians(gq2.type, gq2.ghost_type)(gq2.global_num); Array::const_vector_iterator quad_coords_1 = quadrature_points_coordinates(lq1.type, lq1.ghost_type).begin(spatial_dimension); Array::const_vector_iterator quad_coords_2 = quadrature_points_coordinates(lq2.type, lq2.ghost_type).begin(spatial_dimension); Array & quad_volumes_1 = quadrature_points_volumes(lq1.type, lq1.ghost_type); const Array & jacobians_2 = fem.getIntegratorInterface().getJacobians(gq2.type, gq2.ghost_type); const Real & q2_wJ = jacobians_2(gq2.global_num); // Real & q1_volume = quad_volumes(lq1.global_num); const Vector & q1_coord = quad_coords_1[lq1.global_num]; // quadrature_points_coordinates(lq1.type, lq1.ghost_type).begin(spatial_dimension)[lq1.global_num]; const Vector & q2_coord = quad_coords_2[lq2.global_num]; // quadrature_points_coordinates(lq1.type, lq1.ghost_type).begin(spatial_dimension)[lq2.global_num]; this->weight_func->selectType(lq1.type, lq1.ghost_type, lq2.type, lq2.ghost_type); // Weight function Real r = q1_coord.distance(q2_coord); Real w1 = this->weight_func->operator()(r, lq1, lq2); weight(0) = q2_wJ * w1; // q1_volume += weight(0); quad_volumes_1(lq1.global_num) += weight(0); if(lq2.ghost_type != _ghost && lq1.global_num != lq2.global_num) { const Array & jacobians_1 = fem.getIntegratorInterface().getJacobians(gq1.type, gq1.ghost_type); Array & quad_volumes_2 = quadrature_points_volumes(lq2.type, lq2.ghost_type); const Real & q1_wJ = jacobians_1(gq1.global_num); //Real & q2_volume = quad_volumes_2(lq2.global_num); Real w2 = this->weight_func->operator()(r, lq2, lq1); weight(1) = q1_wJ * w2; quad_volumes_2(lq2.global_num) += weight(1); //q2_volume += weight(1); } else weight(1) = 0.; } } //normalize the weights for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; PairList::const_iterator first_pair = pair_list[ghost_type2].begin(); PairList::const_iterator last_pair = pair_list[ghost_type2].end(); Array::vector_iterator weight_it = pair_weight[ghost_type2]->begin(2); // Compute the weights for(;first_pair != last_pair; ++first_pair, ++weight_it) { Vector & weight = *weight_it; const QuadraturePoint & lq1 = first_pair->first; const QuadraturePoint & lq2 = first_pair->second; Array & quad_volumes_1 = quadrature_points_volumes(lq1.type, lq1.ghost_type); Array & quad_volumes_2 = quadrature_points_volumes(lq2.type, lq2.ghost_type); Real q1_volume = quad_volumes_1(lq1.global_num); weight(0) *= 1. / q1_volume; if(ghost_type2 != _ghost) { Real q2_volume = quad_volumes_2(lq2.global_num); weight(1) *= 1. / q2_volume; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void MaterialNonLocal::weightedAvergageOnNeighbours(const InternalField & to_accumulate, InternalField & accumulated, UInt nb_degree_of_freedom, GhostType ghost_type2) const { AKANTU_DEBUG_IN(); if(ghost_type2 == _not_ghost) { accumulated.reset(); } PairList::const_iterator first_pair = pair_list[ghost_type2].begin(); PairList::const_iterator last_pair = pair_list[ghost_type2].end(); Array::vector_iterator weight_it = pair_weight[ghost_type2]->begin(2); // Compute the weights for(;first_pair != last_pair; ++first_pair, ++weight_it) { Vector & weight = *weight_it; const QuadraturePoint & lq1 = first_pair->first; const QuadraturePoint & lq2 = first_pair->second; const Array & to_acc_1 = to_accumulate(lq1.type, lq1.ghost_type); Array & acc_1 = accumulated(lq1.type, lq1.ghost_type); const Array & to_acc_2 = to_accumulate(lq2.type, lq2.ghost_type); Array & acc_2 = accumulated(lq2.type, lq2.ghost_type); // const Vector & q2_to_acc = to_acc_2[lq2.global_num]; // Vector & q1_acc = acc_1[lq1.global_num]; //q1_acc += weight(0) * q2_to_acc; for(UInt d = 0; d < nb_degree_of_freedom; ++d) { acc_1(lq1.global_num, d) += weight(0) * to_acc_2(lq2.global_num, d); } if(ghost_type2 != _ghost) { for(UInt d = 0; d < nb_degree_of_freedom; ++d) { acc_2(lq2.global_num, d) += weight(1) * to_acc_1(lq1.global_num, d); } } // if(ghost_type2 != _ghost) { // const Vector & q1_to_acc = to_acc_1[lq1.global_num]; // Vector & q2_acc = acc_2[lq2.global_num]; // q2_acc += weight(1) * q1_to_acc; // } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::updateResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); // Update the weights for the non local variable averaging if(ghost_type == _not_ghost && this->weight_func->getUpdateRate() && (this->compute_stress_calls % this->weight_func->getUpdateRate() == 0)) { ElementTypeMapArray quadrature_points_coordinates("quadrature_points_coordinates", getID()); Mesh & mesh = this->model->getFEEngine().getMesh(); mesh.initElementTypeMapArray(quadrature_points_coordinates, spatial_dimension, spatial_dimension); this->fem->computeQuadraturePointsCoordinates(quadrature_points_coordinates, &element_filter); computeWeights(quadrature_points_coordinates); } if(ghost_type == _not_ghost) ++this->compute_stress_calls; computeAllStresses(ghost_type); computeNonLocalStresses(ghost_type); assembleResidual(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::computeAllNonLocalStresses(GhostType ghost_type) { // Update the weights for the non local variable averaging if(ghost_type == _not_ghost) { if(this->weight_func->getUpdateRate() && (this->compute_stress_calls % this->weight_func->getUpdateRate() == 0)) { this->model->getSynchronizerRegistry().asynchronousSynchronize(_gst_mnl_weight); ElementTypeMapArray quadrature_points_coordinates("quadrature_points_coordinates", getID()); Mesh & mesh = this->model->getFEEngine().getMesh(); mesh.initElementTypeMapArray(quadrature_points_coordinates, spatial_dimension, spatial_dimension); this->fem->computeQuadraturePointsCoordinates(quadrature_points_coordinates, &element_filter); this->model->getSynchronizerRegistry().waitEndSynchronize(_gst_mnl_weight); computeWeights(quadrature_points_coordinates); } typename std::map::iterator it = non_local_variables.begin(); typename std::map::iterator end = non_local_variables.end(); for(;it != end; ++it) { NonLocalVariable & non_local_variable = it->second; non_local_variable.non_local->resize(); this->weightedAvergageOnNeighbours(*non_local_variable.local, *non_local_variable.non_local, non_local_variable.nb_component, _not_ghost); } ++this->compute_stress_calls; } else { typename std::map::iterator it = non_local_variables.begin(); typename std::map::iterator end = non_local_variables.end(); for(;it != end; ++it) { NonLocalVariable & non_local_variable = it->second; this->weightedAvergageOnNeighbours(*non_local_variable.local, *non_local_variable.non_local, non_local_variable.nb_component, _ghost); } computeNonLocalStresses(_not_ghost); } } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::savePairs(const std::string & filename) const { std::ofstream pout; std::stringstream sstr; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int prank = comm.whoAmI(); sstr << filename << "." << prank; pout.open(sstr.str().c_str()); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; PairList::const_iterator first_pair = pair_list[ghost_type2].begin(); PairList::const_iterator last_pair = pair_list[ghost_type2].end(); Array::vector_iterator weight_it = pair_weight[ghost_type2]->begin(2); for(;first_pair != last_pair; ++first_pair, ++weight_it) { Vector & weight = *weight_it; const QuadraturePoint & lq1 = first_pair->first; const QuadraturePoint & lq2 = first_pair->second; pout << lq1 << " " << lq2 << " " << weight << std::endl; } } } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::saveWeights(const std::string & filename) const { std::ofstream pout; std::stringstream sstr; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int prank = comm.whoAmI(); sstr << filename << "." << prank; pout.open(sstr.str().c_str()); for (UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type = (GhostType) gt; AKANTU_DEBUG_ASSERT((pair_weight[ghost_type]), "the weights have not been computed yet"); Array & weights = *(pair_weight[ghost_type]); Array::const_vector_iterator weights_it = weights.begin(2); for (UInt i = 0; i < weights.getSize(); ++i, ++weights_it) pout << "w1: " << (*weights_it)(0) <<" w2: " << (*weights_it)(1) << std::endl; } } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::neighbourhoodStatistics(const std::string & filename) const { // std::ofstream pout; // pout.open(filename.c_str()); // const Mesh & mesh = this->model->getFEEngine().getMesh(); // GhostType ghost_type1; // ghost_type1 = _not_ghost; // StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); // Int prank = comm.whoAmI(); // InternalField nb_neighbors("nb_neighbours", *const_cast(this)); // nb_neighbors.initialize(1); // for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { // GhostType ghost_type2 = (GhostType) gt; // UInt existing_pairs_num = gt - _not_ghost; // pair_type::const_iterator first_pair_types = existing_pairs[existing_pairs_num].begin(); // pair_type::const_iterator last_pair_types = existing_pairs[existing_pairs_num].end(); // for (; first_pair_types != last_pair_types; ++first_pair_types) { // const Array & pairs = // pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); // if(prank == 0) { // pout << ghost_type2 << " "; // pout << "Types : " << first_pair_types->first << " " << first_pair_types->second << std::endl; // } // Array::const_iterator< Vector > first_pair = pairs.begin(2); // Array::const_iterator< Vector > last_pair = pairs.end(2); // Array & nb_neigh_1 = nb_neighbors(first_pair_types->first, ghost_type1); // Array & 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 sum_nb_neig = 0; // UInt max_nb_neig = 0; // UInt min_nb_neig = std::numeric_limits::max(); // for(; it != last_type; ++it) { // Array & nb_neighor = nb_neighbors(*it, ghost_type1); // Array::iterator nb_neigh = nb_neighor.begin(); // Array::iterator end_neigh = nb_neighor.end(); // for (; nb_neigh != end_neigh; ++nb_neigh, ++nb_quads) { // UInt nb = *nb_neigh; // sum_nb_neig += nb; // max_nb_neig = std::max(max_nb_neig, nb); // min_nb_neig = std::min(min_nb_neig, nb); // } // } // comm.allReduce(&nb_quads, 1, _so_sum); // comm.allReduce(&sum_nb_neig, 1, _so_sum); // comm.allReduce(&max_nb_neig, 1, _so_max); // comm.allReduce(&min_nb_neig, 1, _so_min); // if(prank == 0) { // pout << ghost_type2 << " "; // pout << "Nb quadrature points: " << nb_quads << std::endl; // Real mean_nb_neig = sum_nb_neig / Real(nb_quads); // pout << ghost_type2 << " "; // pout << "Average nb neighbors: " << mean_nb_neig << "(" << sum_nb_neig << ")" << std::endl; // pout << ghost_type2 << " "; // pout << "Max nb neighbors: " << max_nb_neig << std::endl; // pout << ghost_type2 << " "; // pout << "Min nb neighbors: " << min_nb_neig << std::endl; // } // } // pout.close(); } /* -------------------------------------------------------------------------- */ template inline UInt MaterialNonLocal::getNbDataForElements(const Array & elements, SynchronizationTag tag) const { UInt nb_quadrature_points = this->getModel().getNbQuadraturePoints(elements); UInt size = 0; if(tag == _gst_mnl_for_average) { typename std::map::const_iterator it = non_local_variables.begin(); typename std::map::const_iterator end = non_local_variables.end(); for(;it != end; ++it) { const NonLocalVariable & non_local_variable = it->second; size += non_local_variable.nb_component * sizeof(Real) * nb_quadrature_points; } } size += weight_func->getNbDataForElements(elements, tag); return size; } /* -------------------------------------------------------------------------- */ template inline void MaterialNonLocal::packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const { if(tag == _gst_mnl_for_average) { typename std::map::const_iterator it = non_local_variables.begin(); typename std::map::const_iterator end = non_local_variables.end(); for(;it != end; ++it) { const NonLocalVariable & non_local_variable = it->second; this->packElementDataHelper(*non_local_variable.local, buffer, elements); } } weight_func->packElementData(buffer, elements, tag); } /* -------------------------------------------------------------------------- */ template inline void MaterialNonLocal::unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) { if(tag == _gst_mnl_for_average) { typename std::map::iterator it = non_local_variables.begin(); typename std::map::iterator end = non_local_variables.end(); for(;it != end; ++it) { NonLocalVariable & non_local_variable = it->second; this->unpackElementDataHelper(*non_local_variable.local, buffer, elements); } } weight_func->unpackElementData(buffer, elements, tag); } /* -------------------------------------------------------------------------- */ template inline void MaterialNonLocal::onElementsAdded(const Array & element_list) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ERROR("This is a case not taken into account!!!"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void MaterialNonLocal::onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, __attribute__((unused)) const RemovedElementsEvent & event) { AKANTU_DEBUG_IN(); // Change the pairs in new global numbering for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; PairList::iterator first_pair = pair_list[ghost_type2].begin(); PairList::iterator last_pair = pair_list[ghost_type2].end(); // Array::vector_iterator weight_it = pair_weight[ghost_type2]->begin(2); for(;first_pair != last_pair; ++first_pair) { QuadraturePoint & q1 = first_pair->first; QuadraturePoint gq1 = this->convertToGlobalPoint(q1); q1 = gq1; if(new_numbering.exists(q1.type, q1.ghost_type)) { UInt q1_new_el = new_numbering(q1.type, q1.ghost_type)(gq1.element); AKANTU_DEBUG_ASSERT(q1_new_el != UInt(-1), "A local quadrature_point as been removed instead of just being renumbered"); q1.element = q1_new_el; } QuadraturePoint & q2 = first_pair->second; QuadraturePoint gq2 = this->convertToGlobalPoint(q2); q2 = gq2; if(new_numbering.exists(q2.type, q2.ghost_type)) { UInt q2_new_el = new_numbering(q2.type, q2.ghost_type)(gq2.element); AKANTU_DEBUG_ASSERT(q2_new_el != UInt(-1), "A local quadrature_point as been removed instead of just being renumbered"); q2.element = q2_new_el; } } } // Change the material numbering Material::onElementsRemoved(element_list, new_numbering, event); // Change back the pairs to the new material numbering for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; PairList::iterator first_pair = pair_list[ghost_type2].begin(); PairList::iterator last_pair = pair_list[ghost_type2].end(); // Array::vector_iterator weight_it = pair_weight[ghost_type2]->begin(2); for(;first_pair != last_pair; ++first_pair) { first_pair->first = this->convertToLocalPoint(first_pair->first ); first_pair->second = this->convertToLocalPoint(first_pair->second); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialNonLocal::insertQuadsInNeighborhoods(GhostType ghost_type) { - Real radius = this->weight_func->getRadius(); - const ID & type = this->weight_func->getType(); NonLocalManager & manager = this->model->getNonLocalManager(); UInt spatial_dimension = this->model->getSpatialDimension(); InternalField quadrature_points_coordinates("quadrature_points_coordinates_tmp_nl", *this); quadrature_points_coordinates.initialize(spatial_dimension); /// intialize quadrature point object QuadraturePoint q; q.ghost_type = ghost_type; q.kind = _ek_regular; Mesh::type_iterator it = this->element_filter.firstType(spatial_dimension, ghost_type, _ek_regular); Mesh::type_iterator last_type = this->element_filter.lastType(spatial_dimension, ghost_type, _ek_regular); for(; it != last_type; ++it) { q.type = *it; const Array & elem_filter = this->element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if(nb_element) { UInt nb_quad = this->fem->getNbQuadraturePoints(*it, ghost_type); UInt nb_tot_quad = nb_quad * nb_element; Array & quads = quadrature_points_coordinates(*it, ghost_type); quads.resize(nb_tot_quad); this->model->getFEEngine().computeQuadraturePointsCoordinates(quads, *it, ghost_type, elem_filter); Array::const_vector_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.global_num = q.element * nb_quad + nq; - manager.insertQuad(q, *quad, radius, type); + manager.insertQuad(q, *quad, this->name); ++quad; } ++elem; } } } } + +/* -------------------------------------------------------------------------- */ +template +void MaterialNonLocal::updateNonLocalInternals(ElementTypeMapReal & non_local_flattened, const ID & field_id, const UInt nb_component) { + + for (ghost_type_t::iterator g = ghost_type_t::begin(); g != ghost_type_t::end(); ++g) { + GhostType ghost_type = *g; + /// loop over all types in the material + typedef ElementTypeMapArray:: type_iterator iterator; + iterator it = this->element_filter.firstType(spatial_dimension, ghost_type, _ek_regular); + iterator last_type = this->element_filter.lastType(spatial_dimension, ghost_type, _ek_regular); + for(; it != last_type; ++it) { + ElementType el_type = *it; + Array & internal = this->getInternal(field_id)(el_type, ghost_type); + Array::vector_iterator internal_it = internal.begin(nb_component); + Array & internal_flat = non_local_flattened(el_type, ghost_type); + Array::const_vector_iterator internal_flat_it = internal_flat.begin(nb_component); + /// loop all elements for the given type + const Array & filter = this->element_filter(el_type,ghost_type); + UInt nb_elements = filter.getSize(); + UInt nb_quads = this->getFEEngine().getNbQuadraturePoints(el_type, ghost_type); + for (UInt e = 0; e < nb_elements; ++e) { + UInt global_el = filter(e); + for (UInt q = 0; q < nb_quads; ++q, ++internal_it) { + UInt global_quad = global_el * nb_quads + q; + *internal_it = internal_flat_it[global_quad]; + } + } + } + } +} diff --git a/test/test_model/test_non_local_toolbox/material.dat b/test/test_model/test_non_local_toolbox/material.dat index 6a34fe150..63de6fee8 100644 --- a/test/test_model/test_non_local_toolbox/material.dat +++ b/test/test_model/test_non_local_toolbox/material.dat @@ -1,19 +1,28 @@ material test_material base_wf [ name = mat_1 rho = 1 # density E = 1 # young's modulus nu = 0 # poisson's ratio non_local weight_function [ radius = 0.5 ] ] material test_material base_wf [ name = mat_2 rho = 1 # density E = 1 # young's modulus nu = 0 # poisson's ratio non_local weight_function [ radius = 0.5 ] -] \ No newline at end of file +] + +neighborhoods regions [ + non_local test_region base_wf [ + radius = 0.5 + weight_function weight_parameter [ + update_rate = 1 + ] + ] +] \ No newline at end of file diff --git a/test/test_model/test_non_local_toolbox/material_weight_computation.dat b/test/test_model/test_non_local_toolbox/material_weight_computation.dat index ccef6d2d5..10ca664b5 100644 --- a/test/test_model/test_non_local_toolbox/material_weight_computation.dat +++ b/test/test_model/test_non_local_toolbox/material_weight_computation.dat @@ -1,9 +1,18 @@ material test_material base_wf [ name = mat_1 rho = 1 # density E = 1 # young's modulus nu = 0 # poisson's ratio non_local weight_function [ radius = 0.5 ] ] + +neighborhoods regions [ + non_local test_region base_wf [ + radius = 0.5 + weight_function weight_parameter [ + update_rate = 1 + ] + ] +] \ No newline at end of file diff --git a/test/test_model/test_non_local_toolbox/test_material.cc b/test/test_model/test_non_local_toolbox/test_material.cc index 01171a13a..de846337f 100644 --- a/test/test_model/test_non_local_toolbox/test_material.cc +++ b/test/test_model/test_non_local_toolbox/test_material.cc @@ -1,63 +1,111 @@ /** * @file test_material.cc * @author Aurelia Isabel Cuba Ramos * @date Wed Sep 23 17:16:30 2015 * * @brief Implementation of test material for the non-local neighborhood base test * * @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 "test_material.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template TestMaterial::TestMaterial(SolidMechanicsModel & model, const ID & id) : Material(model, id), MyElasticParent(model, id), MyNonLocalParent(model, id), grad_u_nl("grad_u non local", *this) { this->is_non_local = true; this->grad_u_nl.initialize(dim*dim); } /* -------------------------------------------------------------------------- */ template void TestMaterial::initMaterial() { AKANTU_DEBUG_IN(); this->registerNonLocalVariable(this->gradu, grad_u_nl, spatial_dimension*spatial_dimension); this->model->getNonLocalManager().registerNonLocalVariable(this->gradu.getName(), grad_u_nl.getName(), spatial_dimension*spatial_dimension); MyElasticParent::initMaterial(); MyNonLocalParent::initMaterial(); AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +template +void TestMaterial::insertQuadsInNeighborhoods(GhostType ghost_type) { + + /// this function will add all the quadrature points to the same + /// default neighborhood instead of using one neighborhood per + /// material + NonLocalManager & manager = this->model->getNonLocalManager(); + InternalField quadrature_points_coordinates("quadrature_points_coordinates_tmp_nl", *this); + quadrature_points_coordinates.initialize(spatial_dimension); + + /// intialize quadrature point object + QuadraturePoint q; + q.ghost_type = ghost_type; + q.kind = _ek_regular; + + Mesh::type_iterator it = this->element_filter.firstType(spatial_dimension, ghost_type, _ek_regular); + Mesh::type_iterator last_type = this->element_filter.lastType(spatial_dimension, ghost_type, _ek_regular); + for(; it != last_type; ++it) { + q.type = *it; + const Array & elem_filter = this->element_filter(*it, ghost_type); + UInt nb_element = elem_filter.getSize(); + if(nb_element) { + UInt nb_quad = this->fem->getNbQuadraturePoints(*it, ghost_type); + UInt nb_tot_quad = nb_quad * nb_element; + + Array & quads = quadrature_points_coordinates(*it, ghost_type); + quads.resize(nb_tot_quad); + + this->model->getFEEngine().computeQuadraturePointsCoordinates(quads, *it, ghost_type, elem_filter); + + Array::const_vector_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.global_num = q.element * nb_quad + nq; + manager.insertQuad(q, *quad, "test_region"); + ++quad; + } + ++elem; + } + } + } +} + /* -------------------------------------------------------------------------- */ // Instantiate the material for the 3 dimensions INSTANTIATE_MATERIAL(TestMaterial); /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/test/test_model/test_non_local_toolbox/test_material.hh b/test/test_model/test_non_local_toolbox/test_material.hh index 0bc6ddb80..c7324f37f 100644 --- a/test/test_model/test_non_local_toolbox/test_material.hh +++ b/test/test_model/test_non_local_toolbox/test_material.hh @@ -1,69 +1,71 @@ /** * @file test_material.hh * @author Aurelia Isabel Cuba Ramos * @date Wed Sep 23 17:16:30 2015 * * @brief test material for the non-local neighborhood base test * * @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 "material_elastic.hh" #include "material_non_local.hh" #ifndef __TEST_MATERIAL_HH__ #define __TEST_MATERIAL_HH__ __BEGIN_AKANTU__ template class TestMaterial : public MaterialElastic, public MaterialNonLocal{ /* -------------------------------------------------------------------------- */ /* Constructor/Destructor */ /* -------------------------------------------------------------------------- */ public: TestMaterial(SolidMechanicsModel & model, const ID & id); virtual ~TestMaterial() {}; typedef MaterialNonLocal MyNonLocalParent; typedef MaterialElastic MyElasticParent; /* -------------------------------------------------------------------------- */ /* Methods */ /* -------------------------------------------------------------------------- */ public: void initMaterial(); void computeNonLocalStresses(GhostType ghost_type) {}; + void insertQuadsInNeighborhoods(GhostType ghost_type); + /* -------------------------------------------------------------------------- */ /* Members */ /* -------------------------------------------------------------------------- */ private: InternalField grad_u_nl; }; __END_AKANTU__ #endif /* __TEST_MATERIAL_HH__ */ diff --git a/test/test_model/test_non_local_toolbox/test_non_local_averaging.cc b/test/test_model/test_non_local_toolbox/test_non_local_averaging.cc index fc94f2b03..3d4ec1e53 100644 --- a/test/test_model/test_non_local_toolbox/test_non_local_averaging.cc +++ b/test/test_model/test_non_local_toolbox/test_non_local_averaging.cc @@ -1,76 +1,111 @@ /** * @file test_weight_computation.cc * @author Aurelia Isabel Cuba Ramos * @date Wed Sep 23 16:30:12 2015 * * @brief test for non-local averaging of strain * * @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 "solid_mechanics_model.hh" #include "test_material.hh" #include "non_local_manager.hh" #include "non_local_neighborhood.hh" #include "dumper_paraview.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { - akantu::initialize("material_weight_computation.dat", argc, argv); + akantu::initialize("material.dat", argc, argv); // some configuration variables const UInt spatial_dimension = 2; ElementType element_type = _quadrangle_4; GhostType ghost_type = _not_ghost; // mesh creation and read Mesh mesh(spatial_dimension); mesh.read("plate.msh"); /// model creation SolidMechanicsModel model(mesh); - + + /// creation of material selector + MeshDataMaterialSelector * mat_selector; + mat_selector = new MeshDataMaterialSelector("physical_names", model); + model.setMaterialSelector(*mat_selector); + /// model initialization changed to use our material model.initFull(SolidMechanicsModelOptions(_static, true)); model.registerNewCustomMaterials< TestMaterial >("test_material"); model.initMaterials(); /// dump material index in paraview model.addDumpField("material_index"); + model.addDumpField("grad_u"); + model.addDumpField("grad_u non local"); model.dump(); /// apply constant strain field everywhere in the plate Matrix applied_strain(spatial_dimension, spatial_dimension); applied_strain.clear(); for (UInt i = 0; i < spatial_dimension; ++i) applied_strain(i,i) = 2.; - Array & grad_u = const_cast &>(model.getMaterial(0).getInternal("grad_u")(element_type, ghost_type)); - Array::iterator< Matrix > grad_u_it = grad_u.begin(spatial_dimension, spatial_dimension); - Array::iterator< Matrix > grad_u_end = grad_u.end(spatial_dimension, spatial_dimension); - for (; grad_u_it != grad_u_end; ++grad_u_it) - (*grad_u_it) += applied_strain; + /// apply constant grad_u field in all elements + for (UInt m = 0; m < model.getNbMaterials(); ++m) { + MaterialNonLocal & mat = dynamic_cast & >(model.getMaterial(m)); + + Array & grad_u = const_cast &> (mat.getInternal("grad_u")(element_type, ghost_type)); + Array::iterator< Matrix > grad_u_it = grad_u.begin(spatial_dimension, spatial_dimension); + Array::iterator< Matrix > grad_u_end = grad_u.end(spatial_dimension, spatial_dimension); + for (; grad_u_it != grad_u_end; ++grad_u_it) + (*grad_u_it) += applied_strain; + } /// compute the non-local strains - model.getNonLocalManager().averageInternals(ghost_type); - + model.getNonLocalManager().computeAllNonLocalStresses(); + model.dump(); + + /// verify the result: non-local averaging over constant field must + /// yield same constant field + Real test_result = 0.; + Matrix difference(spatial_dimension, spatial_dimension); + for (UInt m = 0; m < model.getNbMaterials(); ++m) { + MaterialNonLocal & mat = dynamic_cast & >(model.getMaterial(m)); + + Array & grad_u_nl = const_cast &> (mat.getInternal("grad_u non local")(element_type, ghost_type)); + + Array::iterator< Matrix > grad_u_nl_it = grad_u_nl.begin(spatial_dimension, spatial_dimension); + Array::iterator< Matrix > grad_u_nl_end = grad_u_nl.end(spatial_dimension, spatial_dimension); + for (; grad_u_nl_it != grad_u_nl_end; ++grad_u_nl_it) { + difference = (*grad_u_nl_it) - applied_strain; + test_result += difference.norm(); + } + } + + if (test_result > 10.e-13) { + std::cout << "the total norm is: " << test_result << std::endl; + return EXIT_FAILURE; + } + return EXIT_SUCCESS; } diff --git a/test/test_model/test_non_local_toolbox/test_weight_computation.cc b/test/test_model/test_non_local_toolbox/test_weight_computation.cc index ae7b333e3..955e47d9c 100644 --- a/test/test_model/test_non_local_toolbox/test_weight_computation.cc +++ b/test/test_model/test_non_local_toolbox/test_weight_computation.cc @@ -1,72 +1,72 @@ /** * @file test_weight_computation.cc * @author Aurelia Isabel Cuba Ramos * @date Wed Sep 23 16:30:12 2015 * * @brief test for the weight computation with base weight function * * @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 "solid_mechanics_model.hh" #include "test_material.hh" #include "non_local_manager.hh" #include "non_local_neighborhood.hh" #include "dumper_paraview.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { akantu::initialize("material_weight_computation.dat", argc, argv); // some configuration variables const UInt spatial_dimension = 2; // mesh creation and read Mesh mesh(spatial_dimension); mesh.read("plate.msh"); /// model creation SolidMechanicsModel model(mesh); /// model initialization changed to use our material model.initFull(SolidMechanicsModelOptions(_static, true)); model.registerNewCustomMaterials< TestMaterial >("test_material"); model.initMaterials(); /// dump material index in paraview model.addDumpField("material_index"); model.dump(); /// save the weights in a file - NonLocalNeighborhood & neighborhood = dynamic_cast (model.getNonLocalManager().getNeighborhood()); + NonLocalNeighborhood & neighborhood = dynamic_cast &> (model.getNonLocalManager().getNeighborhood()); neighborhood.saveWeights("weights"); /// print results to screen for validation std::ifstream weights; weights.open("weights.0"); std::string current_line; while(getline(weights, current_line)) std::cout << current_line << std::endl; weights.close(); finalize(); return EXIT_SUCCESS; }