diff --git a/src/common/component_libmultiscale.cc b/src/common/component_libmultiscale.cc index b90071f..ff32554 100644 --- a/src/common/component_libmultiscale.cc +++ b/src/common/component_libmultiscale.cc @@ -1,368 +1,372 @@ /* -------------------------------------------------------------------------- */ #include "component_libmultiscale.hh" #include "comm_group.hh" #include "container_array.hh" #include "lm_communicator.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ #define secure(...) secure_call([&]() { __VA_ARGS__; }); /* --------------------------------------------------------------------- */ // OutputContainer /* --------------------------------------------------------------------- */ std::ostream &operator<<(std::ostream &os, OutputContainer &out) { out.get().printself(os); return os; } OutputContainer &OutputContainer::operator=(Component &v) { AutoDispatch::ArgumentAny::operator=(v); this->component = &v; return *this; } UInt OutputContainer::evalRelease() const { UInt release = this->get().getRelease(); if (component) return hash_releases(release, component->evalRelease()); return release; } /* ---------------------------------------------------------------------- */ UInt OutputContainer::getRelease() const { UInt release = this->get().getRelease(); if (component) return hash_releases(release, component->getRelease()); return release; } /* --------------------------------------------------------------------- */ // InputContainer /* --------------------------------------------------------------------- */ InputContainer &InputContainer::operator=(OutputContainer &v) { output_reference = &v; return *this; } /* --------------------------------------------------------------------- */ InputContainer &InputContainer::operator=(OutputContainer &&v) { output_reference = std::make_shared(std::forward(v)); return *this; } OutputContainer &InputContainer::get() const { if (not output_reference.has_value()) throw Component::UnconnectedInput{"input is not connected"}; try { auto *ptr = std::any_cast(output_reference); return *ptr; } catch (std::bad_any_cast &e) { } auto ptr = std::any_cast>(output_reference); return *ptr; } /* --------------------------------------------------------------------- */ OutputContainer &InputContainer::eval() const { OutputContainer &out = this->get(); if (out.component != nullptr) { out.component->compute(); } return out; } /* --------------------------------------------------------------------- */ UInt InputContainer::evalRelease() const { OutputContainer &out = this->get(); return out.evalRelease(); } /* --------------------------------------------------------------------- */ UInt InputContainer::getRelease() const { OutputContainer &out = this->get(); return out.getRelease(); } /* -------------------------------------------------------------------------- */ // Component /* -------------------------------------------------------------------------- */ std::ostream &operator<<(std::ostream &os, Component &comp) { comp.printself(os); return os; } /* --------------------------------------------------------------------- */ Component::Component() : calculated_once(false) {} /* --------------------------------------------------------------------- */ InputContainer & Component::getInternalInput(const std::string &requested_input) { auto &&[input, arg_cont_input] = getContainer(requested_input, inputs); if (not arg_cont_input.has_value()) throw UnconnectedInput{"For component '" + this->getID() + "' input '" + input + "' is not connected"}; return arg_cont_input; } /* --------------------------------------------------------------------- */ OutputContainer &Component::getInput(const std::string &requested_input) { return this->getInternalInput(requested_input).eval(); } /* --------------------------------------------------------------------- */ OutputContainer &Component::getOutput(const std::string &output_name) { auto &&[name, output] = getContainer(output_name, outputs); return output; } /* --------------------------------------------------------------------- */ OutputContainer &Component::evalOutput(const std::string &output_name) { this->compute(); return getOutput(output_name); } /* --------------------------------------------------------------------- */ auto Component::evalOutputs() -> const decltype(outputs) & { this->compute(); return outputs; } /* --------------------------------------------------------------------- */ auto Component::getOutputs() -> const decltype(outputs) & { return outputs; } /* --------------------------------------------------------------------- */ void Component::printself(std::ostream &os) const { for (auto &pair : this->inputs) { os << "input " << pair.first << ":"; if (pair.second.has_value()) os << pair.second.get() << "\n"; else os << "not yet defined/computed\n"; } for (auto &pair : this->outputs) { os << "output " << pair.first << ":"; try { os << pair.second.get() << "\n"; } catch (...) { os << "not yet defined/computed\n"; } } } /* --------------------------------------------------------------------- */ void Component::connect(const std::string &input, Component &comp, const std::string &output) { try { this->connect(input, comp.getOutput(output)); } catch (LibMultiScaleException &e) { if (output == "") inputs[input] = comp; else throw e; } } /* --------------------------------------------------------------------- */ OutputContainer &Component::createOutput(const std::string &output) { outputs[output] = OutputContainer(*this, output); return outputs[output]; } /* --------------------------------------------------------------------- */ InputContainer &Component::createInput(const std::string &input) { inputs[input] = InputContainer(); return inputs[input]; } /* --------------------------------------------------------------------- */ void Component::removeInput(const std::string &input) { inputs.erase(input); } /* --------------------------------------------------------------------- */ void Component::removeOutput(const std::string &output) { outputs.erase(output); } /* --------------------------------------------------------------------- */ void Component::compute() { auto secure_call = [&](auto &&f) { try { f(); + } catch (AutoDispatch::no_argument &e) { + throw Component::UnconnectedInput{e.what()}; } catch (Component::UnconnectedInput &e) { throw e; } catch (Component::NotInCommGroup &e) { throw e; } catch (LibMultiScaleException &e) { LM_FATAL_RE(e, "compute of " << this->getID() << " failed"); } catch (std::exception &e) { LM_FATAL("compute of " << this->getID() << " failed => " << e.what()); } catch (...) { LM_FATAL("compute of " << this->getID() << " failed => unknown exception was raised"); } }; bool need_recompute = this->evalDependencies(); bool amIinGroup = comm_group.amIinGroup(); if (not amIinGroup) { DUMP(this->getID() << ": Do not perform compute as I am not in the right group", DBG_INFO); throw NotInCommGroup{"For component '" + this->getID() + " I am not in the right group"}; calculated_once = true; return; } if (!need_recompute && calculated_once) return; // evalutate the inputs secure(this->evalInputs()); // get the context (comm_group, release) from inputs secure(this->acquireInputsContext()); // call the component code secure(this->compute_make_call()); // sets the context to outputs secure(this->propagateContextToOutputs()); calculated_once = true; DUMP(this->getID() << ": component_called", DBG_MESSAGE); } #undef secure /* -------------------------------------------------------------------------- */ void Component::evalInputs() { for (auto &&[name, inp] : this->inputs) { if (not inp.has_value()) throw Component::UnconnectedInput{"Component(" + this->getID() + "): input(" + name + ") was not created/computed/connected"}; // LM_FATAL("Component(" << this->getID() << "): input(" << name // << ") was not created/computed/connected"); inp.eval(); } } /* -------------------------------------------------------------------------- */ void Component::acquireInputsContext() { if (this->inputs.size() == 0) { this->changeRelease(); } - bool single_comm_groups = true; - // CommGroup *group = nullptr; + std::set comm_groups; for (auto &&[name, inp] : this->inputs) { if (not inp.has_value()) LM_FATAL(name << ": input was not created/computed/connected"); auto &connected_output = inp.get(); if (connected_output.has_argument()) { auto &obj = connected_output.get(); - this->acquireContext(obj); - // single_comm_groups = group == obj.getCommGroup(); + comm_groups.insert(obj.getCommGroup().getID()); } } - if (not single_comm_groups) - this->setCommGroup(Communicator::getGroup("self")); + if (comm_groups.size() > 1) + this->setCommGroup(Communicator::getGroup("all")); + else if (comm_groups.size() == 0) + this->setCommGroup(Communicator::getGroup("self")); + else + this->setCommGroup(Communicator::getGroup(*comm_groups.begin())); UInt release = this->evalRelease(); this->setRelease(release); } /* -------------------------------------------------------------------------- */ void Component::changeRelease() { LMObject::changeRelease(); for (auto &&[name, out] : this->outputs) { if (not out.has_argument()) continue; out->acquireContext(*this); } } /* -------------------------------------------------------------------------- */ bool Component::evalDependencies() { UInt release = this->evalRelease(); if (release != this->getRelease()) return true; return false; } /* -------------------------------------------------------------------------- */ void Component::propagateContextToOutputs() { for (auto &&[name, out] : this->outputs) { if (not out.has_value()) DUMP(name << ": output was not created/computed/connected", DBG_WARNING); try { out.get().acquireContext(*this); } catch (...) { DUMP(name << ": output could not receive context", DBG_WARNING); } } } /* -------------------------------------------------------------------------- */ UInt Component::evalRelease() { if (this->inputs.size() == 0) { return this->getRelease(); } std::vector release_inputs; for (auto &&[name, inp] : this->inputs) { UInt release; if (not inp.has_value()) return LMObject::random_release(); try { release = inp.evalRelease(); } catch (AutoDispatch::no_argument &e) { return LMObject::random_release(); } release_inputs.push_back(release); } UInt release = hash_vector(release_inputs); return release; } /* -------------------------------------------------------------------------- */ void Component::clear() { for (auto &&[name, out] : this->outputs) { this->applyArrayOutput([](auto &output) -> void { output.clear(); }, name); } } /* --------------------------------------------------------------------- */ UInt Component::size(const std::string &name) { try { auto &_output = this->getOutput(name); return _output.size(); } catch (AutoDispatch::no_argument &e) { // do nothing: cannot clean unallocated pointer } catch (AutoDispatch::bad_argument_cast &e) { // do nothing: it is not an array } catch (std::bad_cast &e) { // do nothing: it is not an array } return 0; } /* --------------------------------------------------------------------- */ __END_LIBMULTISCALE__ /* ---------------------------------------------------------------------- */ diff --git a/src/common/lm_object.cc b/src/common/lm_object.cc index fa1f943..d0160ae 100644 --- a/src/common/lm_object.cc +++ b/src/common/lm_object.cc @@ -1,57 +1,56 @@ #include "lm_object.hh" #include "lm_communicator.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ int LMObject::random_release() { // std::srand(std::time(nullptr)); return std::rand(); } /* -------------------------------------------------------------------------- */ LMObject::LMObject(const LMID &id) : id(id), comm_group(Communicator::getGroup("self")) { release = random_release(); } /* -------------------------------------------------------------------------- */ const LMID &LMObject::getID() const { return id; } /* -------------------------------------------------------------------------- */ void LMObject::setID(const LMID &id) { this->id = id; } /* -------------------------------------------------------------------------- */ void LMObject::setCommGroup(CommGroup group) { comm_group = group; } /* -------------------------------------------------------------------------- */ CommGroup LMObject::getCommGroup() const { return comm_group; } /* -------------------------------------------------------------------------- */ void LMObject::acquireContext(const LMObject &obj) { - if (this->comm_group == Communicator::getGroup("self")) - this->setCommGroup(obj.getCommGroup()); + this->setCommGroup(obj.getCommGroup()); void *ptr = this; UInt low = ((UInt *)ptr)[0]; UInt big = ((UInt *)ptr)[1]; this->setRelease(hash_releases(low, big, obj.release)); } /* -------------------------------------------------------------------------- */ UInt LMObject::getRelease() const { return release; } /* -------------------------------------------------------------------------- */ void LMObject::setRelease(UInt r) { release = r; } /* -------------------------------------------------------------------------- */ // void LMObject::incRelease() { ++release; } /* -------------------------------------------------------------------------- */ void LMObject::changeRelease() { release = random_release(); } __END_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ diff --git a/src/coupling/bridging.cc b/src/coupling/bridging.cc index f9f6622..d87f941 100644 --- a/src/coupling/bridging.cc +++ b/src/coupling/bridging.cc @@ -1,416 +1,413 @@ /** * @file bridging.cc * * @author Guillaume Anciaux * * @date Fri Jul 11 15:47:44 2014 * * @brief Bridging object between atomistic and finite elements * * @section LICENSE * * Copyright INRIA and CEA * * The LibMultiScale is a C++ parallel framework for the multiscale * coupling methods dedicated to material simulations. This framework * provides an API which makes it possible to program coupled simulations * and integration of already existing codes. * * This Project was initiated in a collaboration between INRIA Futurs Bordeaux * within ScAlApplix team and CEA/DPTA Ile de France. * The project is now continued at the Ecole Polytechnique Fédérale de Lausanne * within the LSMS/ENAC laboratory. * * This software is governed by the CeCILL-C license under French law and * abiding by the rules of distribution of free software. You can use, * modify and/ or redistribute the software under the terms of the CeCILL-C * license as circulated by CEA, CNRS and INRIA at the following URL * "http://www.cecill.info". * * As a counterpart to the access to the source code and rights to copy, * modify and redistribute granted by the license, users are provided only * with a limited warranty and the software's author, the holder of the * economic rights, and the successive licensors have only limited * liability. * * In this respect, the user's attention is drawn to the risks associated * with loading, using, modifying and/or developing or reproducing the * software by the user in light of its specific status of free software, * that may mean that it is complicated to manipulate, and that also * therefore means that it is reserved for developers and experienced * professionals having in-depth computer knowledge. Users are therefore * encouraged to load and test the software's suitability as regards their * requirements in conditions enabling the security of their systems and/or * data to be ensured and, more generally, to use and operate it in the * same conditions as regards security. * * The fact that you are presently reading this means that you have had * knowledge of the CeCILL-C license and that you accept its terms. * */ /* -------------------------------------------------------------------------- */ #include "bridging.hh" #include "compute_extract.hh" #include "factory_multiscale.hh" #include "filter_geometry.hh" #include "geometry_manager.hh" #include "lm_common.hh" #include "trace_atom.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ Bridging::Bridging(const std::string &name) : LMObject(name), DofAssociation(name), unmatchedPointList("unmatchedPointList-" + name), unmatchedMeshList("unmatchedMeshList-" + name), pointList(createOutput("pointList")), meshList(createOutput("meshList")), contMesh(createOutput("contMesh")), grid(createOutput("grid")) { assoc_found = 0; - unmatchedPointList.setCommGroup(this->comm_group_A); - unmatchedMeshList.setCommGroup(this->comm_group_B); - // contains the elements selected by geometry this->createOutput("unmatchedMeshList") = unmatchedMeshList; // contains the points selected by geometry this->createOutput("unmatchedPointList") = unmatchedPointList; } /* -------------------------------------------------------------------------- */ Bridging::~Bridging() {} /* -------------------------------------------------------------------------- */ void Bridging::filterArray(ContainerArray &array) { UInt index = 0; for (UInt i = 0; i < pointToElement.size(); ++i) if (pointToElement[i] != UINT_MAX) { DUMP("for container, atom " << i << " becomes " << index, DBG_ALL); LM_ASSERT(i < array.size() && index < array.size(), "overflow detected: nmax = " << array.size() << " , index = " << index << " and i = " << i); array(index) = array(i); ++index; } } /* -------------------------------------------------------------------------- */ // void Bridging::cumulPBC(ContainerArray &data) { // UInt npairs [[gnu::unused]] = local_pbc_pairs.size(); // DUMP("Detected " << npairs << " local pairs", DBG_INFO); // for (auto &&pair : local_pbc_pairs) { // UInt i1 = pair.first; // UInt i2 = pair.second; // data(i2) += data(i1); // } // for (auto &&pair : local_pbc_pairs) { // UInt i1 = pair.first; // UInt i2 = pair.second; // data(i1) = data(i2); // } // } /* -------------------------------------------------------------------------- */ void Bridging::mesh2Point(const ContainerArray &data_node, ContainerArray &data_atom) { if (this->local_points == 0) { DUMP("Warning : atomic container is empty:" << " cannot interpolate mesh fields on control points (check " "geometries ?!)", DBG_WARNING); return; } data_atom = smatrix.transpose() * data_node.matrix(); } /* -------------------------------------------------------------------------- */ void Bridging::mesh2Point(FieldType field_type, ContainerArray &data_point) { ComputeExtract extracted_field("ComputeExtract:" + this->getID()); extracted_field.setParam("FIELD", field_type); extracted_field.compute(this->meshList); this->mesh2Point(extracted_field.evalArrayOutput(), data_point); } /* -------------------------------------------------------------------------- */ void Bridging::leastSquarePointData(ContainerArray &, ContainerArray &, UInt) { LM_TOIMPLEMENT; } /* -------------------------------------------------------------------------- */ void Bridging::copySlaveValues(ContainerArray &v) { for (auto &&pair : local_pbc_pairs) { UInt i1 = pair.first; UInt i2 = pair.second; v(i1) = v(i2); } } /* -------------------------------------------------------------------------- */ /* Parallel Methods */ /* -------------------------------------------------------------------------- */ void Bridging::commPoint2Continuum(ContainerArray &buf) { if (this->local_points != 0) { this->distributeVectorA2B("communicate buffer", buf); } } /* -------------------------------------------------------------------------- */ void Bridging::commContinuum2Point(ContainerArray &buf) { if (this->local_points != 0) { this->distributeVectorB2A("communicate buffer", buf); } } /* -------------------------------------------------------------------------- */ void Bridging::synchSumBuffer(ContainerArray &buf) { if (this->local_points != 0) { this->synchronizeVectorBySum("synchsumvector", buf); } } /* -------------------------------------------------------------------------- */ void Bridging::filterRejectedContinuumOwners( std::vector> &unassociated_atoms) { if (this->comm_group_A == this->comm_group_B) return; if (not this->in_group_B()) return; DUMP("local points", DBG_MESSAGE); if (this->local_points == 0) return; // loop through the unassociated atoms declared per each processor // The total list of atom index rejected is constructed UInt offset = 0; UInt cpt_unassociated = 0; std::vector atom_to_remove; for (auto &&[proc, c_w] : enumerate(com_with)) { if (not c_w) continue; DUMP("atomic processor " << proc << " unassociated " << unassociated_atoms[proc].size() << " atoms ", DBG_INFO_STARTUP); for (auto &&unassociated : unassociated_atoms[proc]) { // reconstruct the index for atom declared unassociated UInt unassociated_index = offset + unassociated; atom_to_remove.push_back(unassociated_index); pointToElement[unassociated_index] = UINT_MAX; // auto el_index = this->pointToElement[unassociated_index]; // LM_ASSERT(el_index != UINT_MAX, // "this should not happend !! " << unassociated_index); // // remove the column associated to atom in the shapematrix // this->smatrix.removeAtom(unassociated_index); // // alter the atom-element mapping // pointToElement[unassociated_index] = UINT_MAX; ++cpt_unassociated; DUMP("unassociated atom at position " << unassociated_index << " whereas offset = " << offset << " next offset is " << offset + this->nb_points_per_proc[proc], DBG_INFO_STARTUP); } offset += this->nb_points_per_proc[proc]; } smatrix.removeAtoms(atom_to_remove); // corrects the local_points variable if (cpt_unassociated) { DUMP("removed " << cpt_unassociated << " atoms from unassociation. Now local atoms is " << this->local_points, DBG_MESSAGE); this->local_points -= cpt_unassociated; DUMP("and becomes " << this->local_points, DBG_MESSAGE); } // compute contiguous indexes for atoms (compression of indexes) // the result is a mapping from old indexes to new ones // output: mapping_old_to_new_indexes // std::vector mapping_old_to_new_indexes; // mapping_old_to_new_indexes.assign(pointToElement.size(), UINT_MAX); // UInt indirection = 0; // for (auto &&[at_index, el_index] : enumerate(pointToElement)) { // if (el_index == UINT_MAX) // continue; // mapping_old_to_new_indexes[at_index] = indirection; // ++indirection; // } // LM_ASSERT(indirection == this->local_points, "this should not happend " // << indirection << " " // << this->local_points); // checks that number of mapped atoms equals local atoms variable UInt cpt_tmp = 0; for (auto &&el_index : pointToElement) { if (el_index != UINT_MAX) cpt_tmp++; } if (cpt_tmp != this->local_points) LM_FATAL("biip this should not happend " << cpt_tmp << " " << this->local_points); // // now renumber indirection in shapematrix // indirection = 0; // for (auto &&[at_index, el_index] : enumerate(pointToElement)) { // if (el_index == UINT_MAX) // continue; // UInt new_index = mapping_old_to_new_indexes[at_index]; // LM_ASSERT(new_index != UINT_MAX, "This should not happend"); // smatrix.changeAtomIndirection(at_index, new_index, indirection); // ++indirection; // } // this->smatrix.swapAtomIndirections(); // now i set the duo vector this->createDuoVectorB("FE", pointToElement); } /* -------------------------------------------------------------------------- */ void Bridging::setPBCPairs(std::vector> &pairs) { if (this->is_in_continuum()) pbc_pairs = pairs; } /* -------------------------------------------------------------------------- */ UInt Bridging::getNumberPoints() { LM_TOIMPLEMENT; } /* -------------------------------------------------------------------------- */ UInt Bridging::getNumberElems() { LM_TOIMPLEMENT; } /* -------------------------------------------------------------------------- */ UInt Bridging::getNumberNodes() { LM_TOIMPLEMENT; } /* -------------------------------------------------------------------------- */ void Bridging::attachVector(ContainerArray &tab) { this->pointList.get().attachObject(tab); } /* -------------------------------------------------------------------------- */ void Bridging::updateForMigration() { DUMPBYPROC("update migration for " << this->getID(), DBG_INFO, 0); STARTTIMER("syncMigration resize"); // update local number of atoms in atomic part if (this->in_group_A()) { this->local_points = this->getOutput("pointList").size(); this->positions.resize(this->local_points, spatial_dimension); } STOPTIMER("syncMigration resize"); STARTTIMER("syncMigration duo"); this->synchronizeMigration(this->comm_group_A, this->comm_group_B); STOPTIMER("syncMigration duo"); // pointList.setRelease(this->contA.getRelease()); // meshList.setRelease(this->contB.getRelease()); } /* -------------------------------------------------------------------------- */ /* LMDESC Bridging This class implements a bridging zone where point and elements are associated. For debugging pruposes, several set of DOf are registered to the central system: - unmatched-fe-\${couplerID} : the set of nodes and elements that are contained in the geometry provided by GEOMETRY keyword. - unmatched-point-\${couplerID} : the set of points that are contained in the geometry provided by GEOMETRY keyword. - matched-fe-\${couplerID} : the set of nodes and elements that are containing at least one point of unmatched-point-\${couplerID}. - matched-point-\${couplerID} : the set of points that are contained in at least one element of unmatched-fe-\${couplerID}. */ /* LMHERITANCE filter_geometry dof_association */ void Bridging::declareParams() { this->addSubParsableObject(unmatchedMeshList); DofAssociation::declareParams(); } /* -------------------------------------------------------------------------- */ void Bridging::setPointField(FieldType field, ContainerArray &buffer) { if (not this->is_in_atomic()) return; StimulationField impose_field("impose_point_field"); impose_field.setParam("FIELD", field); impose_field.compute("dofs"_input = this->pointList, "field"_input = buffer); } /* -------------------------------------------------------------------------- */ void Bridging::setPointField(FieldType field, Real val, UInt dim) { ContainerArray data_point(this->local_points, dim); for (auto &&v : data_point.rowwise()) { v = val; } setPointField(field, data_point); } /* -------------------------------------------------------------------------- */ void Bridging::addPointField(FieldType field, ContainerArray &buffer) { StimulationField impose_field("impose_point_field"); impose_field.setParam("FIELD", field); impose_field.setParam("ADDITIVE", true); impose_field.compute("dofs"_input = this->pointList, "field"_input = buffer); } /* -------------------------------------------------------------------------- */ void Bridging::setMeshField(FieldType field, ContainerArray &buffer) { StimulationField impose_field("impose_mesh_field"); impose_field.setParam("FIELD", field); impose_field.compute("dofs"_input = this->meshList, "field"_input = buffer); } /* -------------------------------------------------------------------------- */ void Bridging::addMeshField(FieldType field, ContainerArray &buffer) { StimulationField impose_field("impose_mesh_field"); impose_field.setParam("FIELD", field); impose_field.setParam("ADDITIVE", true); impose_field.compute("dofs"_input = this->meshList, "field"_input = buffer); } /* -------------------------------------------------------------------------- */ void Bridging::setPointField(FieldType field, Eigen::VectorXd val) { ContainerArray data_point(this->local_points, val.rows()); for (auto &&v : data_point.rowwise()) { v = val; } setPointField(field, data_point); } /* -------------------------------------------------------------------------- */ void Bridging::setMeshField(FieldType field, Eigen::VectorXd val) { ContainerArray data_mesh(this->local_points, val.rows()); data_mesh = val; setPointField(field, data_mesh); } /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ diff --git a/src/coupling/bridging_inline_impl.hh b/src/coupling/bridging_inline_impl.hh index 2a5837e..dc480e6 100644 --- a/src/coupling/bridging_inline_impl.hh +++ b/src/coupling/bridging_inline_impl.hh @@ -1,507 +1,511 @@ /** * @file bridging.cc * * @author Guillaume Anciaux * * @date Fri Jul 11 15:47:44 2014 * * @brief Bridging object between atomistic and finite elements * * @section LICENSE * * Copyright INRIA and CEA * * The LibMultiScale is a C++ parallel framework for the multiscale * coupling methods dedicated to material simulations. This framework * provides an API which makes it possible to program coupled simulations * and integration of already existing codes. * * This Project was initiated in a collaboration between INRIA Futurs Bordeaux * within ScAlApplix team and CEA/DPTA Ile de France. * The project is now continued at the Ecole Polytechnique Fédérale de Lausanne * within the LSMS/ENAC laboratory. * * This software is governed by the CeCILL-C license under French law and * abiding by the rules of distribution of free software. You can use, * modify and/ or redistribute the software under the terms of the CeCILL-C * license as circulated by CEA, CNRS and INRIA at the following URL * "http://www.cecill.info". * * As a counterpart to the access to the source code and rights to copy, * modify and redistribute granted by the license, users are provided only * with a limited warranty and the software's author, the holder of the * economic rights, and the successive licensors have only limited * liability. * * In this respect, the user's attention is drawn to the risks associated * with loading, using, modifying and/or developing or reproducing the * software by the user in light of its specific status of free software, * that may mean that it is complicated to manipulate, and that also * therefore means that it is reserved for developers and experienced * professionals having in-depth computer knowledge. Users are therefore * encouraged to load and test the software's suitability as regards their * requirements in conditions enabling the security of their systems and/or * data to be ensured and, more generally, to use and operate it in the * same conditions as regards security. * * The fact that you are presently reading this means that you have had * knowledge of the CeCILL-C license and that you accept its terms. * */ /* -------------------------------------------------------------------------- */ #include "bridging.hh" #include "compute_extract.hh" #include "factory_multiscale.hh" #include "filter_geometry.hh" #include "geometry_manager.hh" #include "lm_common.hh" #include "stimulation_field.hh" #include "trace_atom.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ template void Bridging::init(ContA &contA, ContC &contC) { auto all_group = Communicator::getCommunicator().getGroup("all"); DofAssociation::init(contA, contC); const UInt Dim = spatial_dimension; contMesh = contC; + // unmatchedPointList.setCommGroup(this->comm_group_A); + // unmatchedMeshList.setCommGroup(this->comm_group_B); + unmatchedPointList.acquireContext(contA); + unmatchedMeshList.acquireContext(contC); if (in_group_A()) { pointList.alloc("pointList-" + this->getID()); pointList->acquireContext(contA); unmatchedPointList.acquireContext(contA); } if (in_group_B()) { meshList.alloc("meshList-" + this->getID()); meshList->acquireContext(contC); unmatchedMeshList.acquireContext(contC); } buffer_for_points.acquireContext(contA); buffer_for_nodes.acquireContext(contC); DUMPBYPROC("selecting DOFs in bridging zone", DBG_INFO_STARTUP, 0); if (in_group_A()) { this->buildPointList(contA); this->buildPositions(unmatchedPointList.evalOutput()); } if (in_group_B()) this->buildContinuumDOFsList(contC); DUMP(this->getID() << " : Generating Communication Scheme", DBG_INFO_STARTUP); all_group.synchronize(); DUMP(this->getID() << " : exchange coarse geometries", DBG_MESSAGE); if (this->in_group_A()) { this->exchangeGeometries(this->unmatchedPointList.evalOutput()); } if (in_group_B()) { this->exchangeGeometries(this->unmatchedMeshList.evalOutput()); } all_group.synchronize(); DUMP(this->getID() << " : exchange points positions", DBG_MESSAGE); this->exchangePositions(Dim); all_group.synchronize(); DUMP(this->getID() << " : build shape matrix", DBG_MESSAGE); if (in_group_B()) { DUMP(this->getID() << " : build shape matrix (parsing " << this->local_points << " points by position)", DBG_INFO); this->buildShapeMatrix(this->local_points, this->unmatchedMeshList.evalOutput()); this->local_points = this->assoc_found; DUMP(this->getID() << " : I will manage " << this->local_points << " points", DBG_INFO); } all_group.synchronize(); DUMP(this->getID() << " : ExchangeAssociation", DBG_MESSAGE); ContainerArray managed; this->exchangeAssociationInformation(managed, pointToElement); for (auto &&[i, c_w] : enumerate(this->com_with)) DUMP("second test Do I communicate with " << i << " : " << c_w, DBG_INFO); all_group.synchronize(); DUMP(this->getID() << " : generate duo vector and detect double point assignement", DBG_MESSAGE); std::vector> unassociated; if (this->in_group_A()) this->createDuoVectorA("Point", managed, pointToElement, unassociated); DUMP(this->getID() << " : doing exchange rejected", DBG_MESSAGE); this->exchangeRejectedContinuumOwners(unassociated); all_group.synchronize(); DUMP(this->getID() << " : doing filter rejected", DBG_MESSAGE); if (this->in_group_B()) this->filterRejectedContinuumOwners(unassociated); all_group.synchronize(); DUMP(this->getID() << " : filtering position & association vector", DBG_MESSAGE); if (this->in_group_A()) { filterPointListForUnmatched(pointList); } if (this->in_group_B()) { this->filterArray(this->positions); this->filterMapping(pointToElement); // this->buildNodeShape(this->meshList); } DUMP(this->getID() << " : all done for scheme generation", DBG_MESSAGE); all_group.synchronize(); } /* -------------------------------------------------------------------------- */ template void Bridging::buildContinuumDOFsList(ContC &contMesh) { constexpr UInt Dim = ContC::ContainerSubset::Dim; unmatchedMeshList.setParam("GEOMETRY", this->geomID); unmatchedMeshList.compute(contMesh); auto &unmatched_mesh = unmatchedMeshList.evalOutput(); DUMP("we have found " << unmatched_mesh.getContainerElems().size() << " concerned elements", DBG_INFO_STARTUP); STARTTIMER("Filling spatial-grid"); auto contElems = unmatched_mesh.getContainerElems(); auto contNodes = unmatched_mesh.getContainerNodes(); std::vector nodes; Geometry &geom = *GeometryManager::getManager().getGeometry(this->geomID); Cube cube = geom.getBoundingBox(); DUMP("geometry of the grid => \n" << cube, DBG_INFO_STARTUP); auto &grid = this->grid.alloc>(cube, this->grid_division); UInt nb_elem = 0; for (auto &&el : contElems) { std::vector> node_coords; auto &&nodes = el.globalIndexes(); UInt nb = nodes.size(); for (UInt i = 0; i < nb; ++i) { auto nd = el.getNode(i); #ifndef LM_OPTIMIZED Real mass_node = nd.mass(); LM_ASSERT(mass_node > 0, "invalid mass" << nodes[i] << " " << mass_node); #endif node_coords.push_back(nd.position0()); } grid.addFiniteElement(nb_elem, node_coords); ++nb_elem; } STOPTIMER("Filling spatial-grid"); DUMP("we have found " << contElems.size() << " concerned elements", DBG_INFO_STARTUP); DUMP("in average " << grid.getAverageEltByBlock() << " elements per block", DBG_INFO_STARTUP); } /* -------------------------------------------------------------------------- */ template void Bridging::buildPointList(ContA &contPoints) { unmatchedPointList.setParam("GEOMETRY", this->geomID); unmatchedPointList.compute(contPoints); this->local_points = unmatchedPointList.size(); if (!this->local_points) DUMP("We found no point in the bridging zone", DBG_INFO); DUMP("We found " << this->local_points << " point in the bridging zone", DBG_INFO_STARTUP); } /* -------------------------------------------------------------------------- */ template void Bridging::buildShapeMatrix(UInt nb_points, ContMesh &unmatchedMeshList) { pointToElement.assign(nb_points, UINT_MAX); if (unmatchedMeshList.getContainerElems().size() == 0) { // need to free the output produced by the component DUMP("elem_rec is empty!", DBG_WARNING); this->getOutput("grid").reset(); this->removeOutput("grid"); return; } constexpr UInt Dim = ContMesh::ContainerSubset::Dim; auto &grid = this->grid.get &>(); //**************************************************************** /* building the association array points <-> elements */ //**************************************************************** STARTTIMER("Construction association"); for (UInt i = 0; i < nb_points; ++i) { #ifndef TIMER if (i % 100000 == 0) DUMP("building association: point " << i << "/" << nb_points, DBG_INFO); #endif Vector x = this->positions(i); pointToElement[i] = UINT_MAX; for (auto &&j : grid.findSet(x)) { auto el = unmatchedMeshList.getContainerElems().get(j); DUMP("is point " << i << " within element " << j << "?", DBG_ALL); if (el.contains(x)) { LM_ASSERT(pointToElement[i] == UINT_MAX, "this should not happen. " << "I found twice the same point while filtering"); DUMP("associating point " << i << " and element " << j, DBG_ALL); pointToElement[i] = j; ++assoc_found; break; } } } STOPTIMER("Construction association"); std::vector nb_points_per_element( unmatchedMeshList.getContainerElems().size()); /* build number of atoms per element */ for (UInt i = 0; i < nb_points; ++i) { if (pointToElement[i] != UINT_MAX) ++nb_points_per_element[pointToElement[i]]; } auto &meshList = this->getOutput("meshList"); /* filter element&node list for the unassociated nodes&elements */ filterContainerElems(nb_points_per_element, this->contMesh); auto &elements = meshList.getContainerElems(); auto &nodes = meshList.getContainerNodes(); smatrix.resize(nodes.size(), assoc_found); /* set the indirection (alloc tab) plus the shapes */ std::vector shapes; UInt atom_index = 0; for (UInt i = 0; i < nb_points; ++i) { if (pointToElement[i] == UINT_MAX) continue; auto &&X = this->positions(i); auto el = elements.get(pointToElement[i]); el.computeShapes(shapes, X); auto &&node_indices = el.getAlteredConnectivity(); IF_TRACED(X, "checking for atom in trouble"); std::vector local_nodes(shapes.size()); for (UInt nd = 0; nd < shapes.size(); ++nd) smatrix.insert(node_indices[nd], atom_index) = shapes[nd]; ++atom_index; } } /* -------------------------------------------------------------------------- */ template void Bridging::buildLocalPBCPairs(ContC &meshList) { UInt npairs = pbc_pairs.size(); for (UInt pair = 0; pair < npairs; ++pair) { UInt ind1 = pbc_pairs[pair].first; UInt ind2 = pbc_pairs[pair].second; UInt i1 = meshList.subIndex2Index(ind1); UInt i2 = meshList.subIndex2Index(ind2); if (i1 == UINT_MAX || i2 == UINT_MAX) continue; std::pair p((UInt)(i1), (UInt)(i2)); local_pbc_pairs.push_back(p); } DUMP("Detected " << local_pbc_pairs.size() << " local pairs", DBG_MESSAGE); LM_TOIMPLEMENT; // smatrix->alterMatrixIndexesForPBC(local_pbc_pairs); } /* -------------------------------------------------------------------------- */ // template void Bridging::buildNodeShape(ContC &meshList) { // if (this->local_points == 0) // return; // UInt nb_coupled_nodes = meshList.size(); // node_shape.assign(nb_coupled_nodes, 0); // LM_ASSERT(smatrix, "shapematrix object was not allocated " // << " eventhough local atoms were detected"); // smatrix->buildNodeShape(node_shape); // buildLocalPBCPairs(this->meshList); // cumulPBC(node_shape); // #ifndef LM_OPTIMIZED // for (UInt i = 0; i < nb_coupled_nodes; ++i) { // if (node_shape[i] <= 1e-1) { // DUMP(i << " " << node_shape[i], DBG_MESSAGE); // LM_FATAL("Please check your geometry:" // << " one or more nodes has zero nodal shape values." // << " Changing the FE element size may resolve this issue!" // << std::endl); // } // } // #endif // LM_OPTIMIZED // FilterManager::getManager().addObject(node_shape); // } /* -------------------------------------------------------------------------- */ template void Bridging::filterContainerElems(std::vector &nb_atom_per_element, ContMesh &contMesh) { using ContainerSubset = typename ContMesh::ContainerSubset; typename ContainerSubset::ContainerElems new_t(this->getID() + ":elems"); std::vector new_index; auto &unmatchedMeshList = this->unmatchedMeshList.evalOutput(); auto &meshList = this->meshList.get(); for (UInt el = 0; el < nb_atom_per_element.size(); ++el) { if (nb_atom_per_element[el] > 0) { DUMP("elem " << el << " becomes " << new_t.size(), DBG_ALL); new_index.push_back(new_t.size()); nb_atom_per_element[new_t.size()] = nb_atom_per_element[el]; new_t.push_back(unmatchedMeshList.getContainerElems().get(el)); } else new_index.push_back(new_t.size()); } nb_atom_per_element.resize(new_t.size()); meshList.clear(); // rebuild the element container meshList.getContainerElems() = new_t; // rebuild the association list for (UInt i = 0; i < pointToElement.size(); ++i) { if (pointToElement[i] != UINT_MAX) { pointToElement[i] = new_index[pointToElement[i]]; } } meshList.computeAlteredConnectivity(contMesh.getContainerNodes()); UInt nb = meshList.size(); if (!nb) { DUMP("We found no nodes in bridging zone", DBG_INFO_STARTUP); return; } DUMP("We found " << nb << " nodes concerned into " << meshList.getContainerElems().size() << " elements", DBG_INFO_STARTUP); } /* -------------------------------------------------------------------------- */ template void Bridging::filterPointListForUnmatched(ContA &pointList) { auto &unmatchedPointList = this->unmatchedPointList.evalOutput(); pointList.clear(); for (UInt i = 0; i < pointToElement.size(); ++i) if (pointToElement[i] != UINT_MAX) { DUMP("for container, atom " << i << " becomes " << pointList.size(), DBG_ALL); pointList.push_back(unmatchedPointList.get(i)); } else { DUMP("atom filtered " << i << " " << unmatchedPointList.get(i), DBG_DETAIL); } pointList.acquireContext(unmatchedPointList); } /* -------------------------------------------------------------------------- */ template void Bridging::solveLeastSquare(ContainerArray &mesh_data, ContainerArray &atomic_data, ContC &meshList) { // build right hand side of leastsquare system mesh_data.resize(meshList.size(), atomic_data.cols()); mesh_data = smatrix * atomic_data.matrix(); cumulPBC(mesh_data); LM_TOIMPLEMENT; // // solve least square system // for (UInt i = 0; i < meshList.size(); ++i) { // mesh_data(i) /= node_shape[i]; // } } /* -------------------------------------------------------------------------- */ // template // void Bridging::averageOnElements(ContainerArray &data_atomic, // ContainerArray &data_mesh, // ContC &meshList [[gnu::unused]]) { // LM_TOIMPLEMENT; // // data_mesh.assign(meshList.getContainerElems().rows(), // // meshList.getContainerElems().cols()); // // smatrix->averageOnElements(data_atomic, data_mesh); // } /* -------------------------------------------------------------------------- */ // DECLARE_BRIDGING_ATOMIC_CONTINUUM_TEMPLATE(Bridging) // DECLARE_BRIDGING_DD_CONTINUUM_TEMPLATE(Bridging) // template // void Bridging::buildLeastSquareMatrix(SparseMatrix &mat) { // mat = smatrix * smatrix.transpose(); // cumulPBC(mat); // } // /* -------------------------------------------------------------------------- // */ template void Bridging::cumulPBC(Matrix &mat) { for (auto &&[i1, i2] : local_pbc_pairs) { mat.row(i2) += mat.row(i1); } for (auto &&[i1, i2] : local_pbc_pairs) { mat.row(i1) += mat.row(i2); } } /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__