Diffusion akantu (features/34-poisson-solver-with-multiple-dofs) Browse

34-poisson…

dof_manager_default.cc
View Options

			/**
			* @file dof_manager_default.cc
			*
			* @author Nicolas Richart <nicolas.richart@epfl.ch>
			*
			* @date creation: Tue Aug 18 2015
			* @date last modification: Thu Feb 08 2018
			*
			* @brief Implementation of the default DOFManager
			*
			* @section LICENSE
			*
			* Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
			* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
			*
			* Akantu is free software: you can redistribute it and/or modify it under the
			* terms of the GNU Lesser General Public License as published by the Free
			* Software Foundation, either version 3 of the License, or (at your option) any
			* later version.
			*
			* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
			* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
			* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
			* details.
			*
			* You should have received a copy of the GNU Lesser General Public License
			* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
			*
			*/

			/* -------------------------------------------------------------------------- */
			#include "dof_manager_default.hh"
			#include "communicator.hh"
			#include "dof_synchronizer.hh"
			#include "element_group.hh"
			#include "node_synchronizer.hh"
			#include "non_linear_solver_default.hh"
			#include "periodic_node_synchronizer.hh"
			#include "solver_vector_default.hh"
			#include "solver_vector_distributed.hh"
			#include "sparse_matrix_aij.hh"
			#include "terms_to_assemble.hh"
			#include "time_step_solver_default.hh"
			/* -------------------------------------------------------------------------- */
			#include <algorithm>
			#include <memory>
			#include <numeric>
			#include <unordered_map>
			/* -------------------------------------------------------------------------- */

			namespace akantu {

			/* -------------------------------------------------------------------------- */
			inline void DOFManagerDefault::addSymmetricElementalMatrixToSymmetric(
			SparseMatrixAIJ & matrix, const Matrix<Real> & elementary_mat,
			const Vector<Int> & equation_numbers, UInt max_size) {
			for (UInt i = 0; i < elementary_mat.rows(); ++i) {
			UInt c_irn = equation_numbers(i);
			if (c_irn < max_size) {
			for (UInt j = i; j < elementary_mat.cols(); ++j) {
			UInt c_jcn = equation_numbers(j);
			if (c_jcn < max_size) {
			matrix(c_irn, c_jcn) += elementary_mat(i, j);
			}
			}
			}
			}
			}

			/* -------------------------------------------------------------------------- */
			inline void DOFManagerDefault::addUnsymmetricElementalMatrixToSymmetric(
			SparseMatrixAIJ & matrix, const Matrix<Real> & elementary_mat,
			const Vector<Int> & equation_numbers, UInt max_size) {
			for (UInt i = 0; i < elementary_mat.rows(); ++i) {
			UInt c_irn = equation_numbers(i);
			if (c_irn < max_size) {
			for (UInt j = 0; j < elementary_mat.cols(); ++j) {
			UInt c_jcn = equation_numbers(j);
			if (c_jcn < max_size) {
			if (c_jcn >= c_irn) {
			matrix(c_irn, c_jcn) += elementary_mat(i, j);
			}
			}
			}
			}
			}
			}

			/* -------------------------------------------------------------------------- */
			inline void DOFManagerDefault::addElementalMatrixToUnsymmetric(
			SparseMatrixAIJ & matrix, const Matrix<Real> & elementary_mat,
			const Vector<Int> & equation_numbers, UInt max_size) {
			for (UInt i = 0; i < elementary_mat.rows(); ++i) {
			UInt c_irn = equation_numbers(i);
			if (c_irn < max_size) {
			for (UInt j = 0; j < elementary_mat.cols(); ++j) {
			UInt c_jcn = equation_numbers(j);
			if (c_jcn < max_size) {
			matrix(c_irn, c_jcn) += elementary_mat(i, j);
			}
			}
			}
			}
			}

			/* -------------------------------------------------------------------------- */
			DOFManagerDefault::DOFManagerDefault(const ID & id, const MemoryID & memory_id)
			: DOFManager(id, memory_id),
			global_blocked_dofs(0, 1, std::string(id + ":global_blocked_dofs")),
			previous_global_blocked_dofs(
			0, 1, std::string(id + ":previous_global_blocked_dofs")),
			dofs_flag(0, 1, std::string(id + ":dofs_type")),
			data_cache(0, 1, std::string(id + ":data_cache_array")),
			global_equation_number(0, 1, "global_equation_number"),
			synchronizer(nullptr) {
			residual = std::make_unique<SolverVectorDefault>(
			*this, std::string(id + ":residual"));
			solution = std::make_unique<SolverVectorDefault>(
			*this, std::string(id + ":solution"));
			}

			/* -------------------------------------------------------------------------- */
			DOFManagerDefault::DOFManagerDefault(Mesh & mesh, const ID & id,
			const MemoryID & memory_id)
			: DOFManager(mesh, id, memory_id),
			global_blocked_dofs(0, 1, std::string(id + ":global_blocked_dofs")),
			previous_global_blocked_dofs(
			0, 1, std::string(id + ":previous_global_blocked_dofs")),
			dofs_flag(0, 1, std::string(id + ":dofs_type")),
			data_cache(0, 1, std::string(id + ":data_cache_array")),
			global_equation_number(0, 1, "global_equation_number"),
			first_global_dof_id(0), synchronizer(nullptr) {
			if (this->mesh->isDistributed()) {
			this->synchronizer = std::make_unique<DOFSynchronizer>(
			*this, this->id + ":dof_synchronizer", this->memory_id);
			residual = std::make_unique<SolverVectorDistributed>(
			*this, std::string(id + ":residual"));
			solution = std::make_unique<SolverVectorDistributed>(
			*this, std::string(id + ":solution"));
			} else {
			residual = std::make_unique<SolverVectorDefault>(
			*this, std::string(id + ":residual"));
			solution = std::make_unique<SolverVectorDefault>(
			*this, std::string(id + ":solution"));
			}
			}

			/* -------------------------------------------------------------------------- */
			DOFManagerDefault::~DOFManagerDefault() = default;

			/* -------------------------------------------------------------------------- */
			template <typename T>
			void DOFManagerDefault::makeConsistentForPeriodicity(const ID & dof_id,
			Array<T> & array) {
			auto & dof_data = this->getDOFDataTyped<DOFDataDefault>(dof_id);
			if (dof_data.support_type != _dst_nodal)
			return;

			if (not mesh->isPeriodic())
			return;

			this->mesh->getPeriodicNodeSynchronizer()
			.reduceSynchronizeWithPBCSlaves<AddOperation>(array);
			}

			/* -------------------------------------------------------------------------- */
			template <typename T>
			void DOFManagerDefault::assembleToGlobalArray(
			const ID & dof_id, const Array<T> & array_to_assemble,
			Array<T> & global_array, T scale_factor) {
			AKANTU_DEBUG_IN();

			auto & dof_data = this->getDOFDataTyped<DOFDataDefault>(dof_id);
			AKANTU_DEBUG_ASSERT(dof_data.local_equation_number.size() ==
			array_to_assemble.size() *
			array_to_assemble.getNbComponent(),
			"The array to assemble does not have a correct size."
			<< " (" << array_to_assemble.getID() << ")");

			if (dof_data.support_type == _dst_nodal and mesh->isPeriodic()) {
			for (auto && data :
			zip(dof_data.local_equation_number, dof_data.associated_nodes,
			make_view(array_to_assemble))) {
			auto && equ_num = std::get<0>(data);
			auto && node = std::get<1>(data);
			auto && arr = std::get<2>(data);
			global_array(equ_num) +=
			scale_factor * (arr) * (not this->mesh->isPeriodicSlave(node));
			}
			} else {
			for (auto && data :
			zip(dof_data.local_equation_number, make_view(array_to_assemble))) {
			auto && equ_num = std::get<0>(data);
			auto && arr = std::get<1>(data);
			global_array(equ_num) += scale_factor * (arr);
			}
			}
			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			DOFManagerDefault::DOFDataDefault::DOFDataDefault(const ID & dof_id)
			: DOFData(dof_id), associated_nodes(0, 1, dof_id + "associated_nodes") {}

			/* -------------------------------------------------------------------------- */
			DOFManager::DOFData & DOFManagerDefault::getNewDOFData(const ID & dof_id) {
			this->dofs[dof_id] = std::make_unique<DOFDataDefault>(dof_id);
			return *this->dofs[dof_id];
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::registerDOFsInternal(const ID & dof_id, UInt nb_dofs,
			UInt nb_pure_local_dofs) {
			// access the relevant data to update
			auto & dof_data = this->getDOFDataTyped<DOFDataDefault>(dof_id);
			const auto & support_type = dof_data.support_type;
			const auto & group = dof_data.group_support;

			switch (support_type) {
			case _dst_nodal:
			if (group != "__mesh__") {
			auto & support_nodes =
			this->mesh->getElementGroup(group).getNodeGroup().getNodes();
			this->updateDOFsData(
			dof_data, nb_dofs, nb_pure_local_dofs, support_nodes.size(),
			[&support_nodes](UInt node) -> UInt { return support_nodes[node]; });
			} else {
			this->updateDOFsData(dof_data, nb_dofs, nb_pure_local_dofs,
			mesh->getNbNodes(),
			[](UInt node) -> UInt { return node; });
			}
			break;
			case _dst_generic:
			this->updateDOFsData(dof_data, nb_dofs, nb_pure_local_dofs);
			break;
			}

			// update the synchronizer if needed
			if (this->synchronizer)
			this->synchronizer->registerDOFs(dof_id);
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::registerDOFs(const ID & dof_id,
			Array<Real> & dofs_array,
			const DOFSupportType & support_type) {
			// stores the current numbers of dofs
			UInt nb_dofs_old = this->local_system_size;
			UInt nb_pure_local_dofs_old = this->pure_local_system_size;

			// update or create the dof_data
			DOFManager::registerDOFs(dof_id, dofs_array, support_type);

			UInt nb_dofs = this->local_system_size - nb_dofs_old;
			UInt nb_pure_local_dofs =
			this->pure_local_system_size - nb_pure_local_dofs_old;

			this->registerDOFsInternal(dof_id, nb_dofs, nb_pure_local_dofs);
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::registerDOFs(const ID & dof_id,
			Array<Real> & dofs_array,
			const ID & group_support) {
			// stores the current numbers of dofs
			UInt nb_dofs_old = this->local_system_size;
			UInt nb_pure_local_dofs_old = this->pure_local_system_size;

			// update or create the dof_data
			DOFManager::registerDOFs(dof_id, dofs_array, group_support);

			UInt nb_dofs = this->local_system_size - nb_dofs_old;
			UInt nb_pure_local_dofs =
			this->pure_local_system_size - nb_pure_local_dofs_old;

			this->registerDOFsInternal(dof_id, nb_dofs, nb_pure_local_dofs);
			}

			/* -------------------------------------------------------------------------- */
			SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & id,
			const MatrixType & matrix_type) {
			return this->registerSparseMatrix<SparseMatrixAIJ>(*this, id, matrix_type);
			}

			/* -------------------------------------------------------------------------- */
			SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & id,
			const ID & matrix_to_copy_id) {
			return this->registerSparseMatrix<SparseMatrixAIJ>(id, matrix_to_copy_id);
			}

			/* -------------------------------------------------------------------------- */
			SparseMatrixAIJ & DOFManagerDefault::getMatrix(const ID & id) {
			auto & matrix = DOFManager::getMatrix(id);
			return dynamic_cast<SparseMatrixAIJ &>(matrix);
			}

			/* -------------------------------------------------------------------------- */
			NonLinearSolver &
			DOFManagerDefault::getNewNonLinearSolver(const ID & id,
			const NonLinearSolverType & type) {
			switch (type) {
			#if defined(AKANTU_IMPLICIT)
			case NonLinearSolverType::_newton_raphson:
			/* FALLTHRU */
			/* [[fallthrough]]; un-comment when compiler will get it */
			case NonLinearSolverType::_newton_raphson_modified: {
			return this->registerNonLinearSolver<NonLinearSolverNewtonRaphson>(
			*this, id, type);
			}
			case NonLinearSolverType::_linear: {
			return this->registerNonLinearSolver<NonLinearSolverLinear>(*this, id,
			type);
			}
			#endif
			case NonLinearSolverType::_lumped: {
			return this->registerNonLinearSolver<NonLinearSolverLumped>(*this, id,
			type);
			}
			default:
			AKANTU_EXCEPTION("The asked type of non linear solver is not supported by "
			"this dof manager");
			}
			}

			/* -------------------------------------------------------------------------- */
			TimeStepSolver &
			DOFManagerDefault::getNewTimeStepSolver(const ID & id,
			const TimeStepSolverType & type,
			NonLinearSolver & non_linear_solver) {
			return this->registerTimeStepSolver<TimeStepSolverDefault>(*this, id, type,
			non_linear_solver);
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::clearResidual() {
			this->residual->resize();
			this->residual->clear();
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::clearMatrix(const ID & mtx) {
			this->getMatrix(mtx).clear();
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::clearLumpedMatrix(const ID & mtx) {
			this->getLumpedMatrix(mtx).clear();
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::updateGlobalBlockedDofs() {
			auto it = this->dofs.begin();
			auto end = this->dofs.end();

			this->previous_global_blocked_dofs.copy(this->global_blocked_dofs);
			this->global_blocked_dofs.resize(this->local_system_size);
			this->global_blocked_dofs.clear();

			for (; it != end; ++it) {
			if (!this->hasBlockedDOFs(it->first))
			continue;

			DOFData & dof_data = *it->second;
			this->assembleToGlobalArray(it->first, *dof_data.blocked_dofs,
			this->global_blocked_dofs, true);
			}
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::splitSolutionPerDOFs() {
			auto it = this->dofs.begin();
			auto end = this->dofs.end();

			for (; it != end; ++it) {
			DOFData & dof_data = *it->second;
			dof_data.solution.resize(dof_data.dof->size() *
			dof_data.dof->getNbComponent());
			this->getSolutionPerDOFs(it->first, dof_data.solution);
			}
			}

			/* -------------------------------------------------------------------------- */
			template <typename T>
			void DOFManagerDefault::getArrayPerDOFs(const ID & dof_id,
			const Array<T> & global_array,
			Array<T> & local_array) const {
			AKANTU_DEBUG_IN();

			const Array<Int> & equation_number = this->getLocalEquationsNumbers(dof_id);

			UInt nb_degree_of_freedoms = equation_number.size();
			local_array.resize(nb_degree_of_freedoms / local_array.getNbComponent());

			auto loc_it = local_array.begin_reinterpret(nb_degree_of_freedoms);
			auto equ_it = equation_number.begin();

			for (UInt d = 0; d < nb_degree_of_freedoms; ++d, ++loc_it, ++equ_it) {
			(loc_it) = global_array(equ_it);
			}

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::getSolutionPerDOFs(const ID & dof_id,
			Array<Real> & solution_array) {
			AKANTU_DEBUG_IN();
			this->getArrayPerDOFs(dof_id, this->getSolutionArray(), solution_array);
			AKANTU_DEBUG_OUT();
			}
			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::getLumpedMatrixPerDOFs(const ID & dof_id,
			const ID & lumped_mtx,
			Array<Real> & lumped) {
			AKANTU_DEBUG_IN();
			this->getArrayPerDOFs(dof_id, this->getLumpedMatrix(lumped_mtx), lumped);
			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::assembleToResidual(const ID & dof_id,
			Array<Real> & array_to_assemble,
			Real scale_factor) {
			AKANTU_DEBUG_IN();

			this->makeConsistentForPeriodicity(dof_id, array_to_assemble);

			this->assembleToGlobalArray(dof_id, array_to_assemble,
			this->getResidualArray(), scale_factor);

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::assembleToLumpedMatrix(const ID & dof_id,
			Array<Real> & array_to_assemble,
			const ID & lumped_mtx,
			Real scale_factor) {
			AKANTU_DEBUG_IN();

			this->makeConsistentForPeriodicity(dof_id, array_to_assemble);

			Array<Real> & lumped = this->getLumpedMatrix(lumped_mtx);
			this->assembleToGlobalArray(dof_id, array_to_assemble, lumped, scale_factor);

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::assembleMatMulVectToGlobalArray(const ID & dof_id,
			const ID & A_id,
			const Array<Real> & x,
			Array<Real> & array,
			Real scale_factor) {
			SparseMatrixAIJ & A = this->getMatrix(A_id);

			// Array<Real> data_cache(this->local_system_size, 1, 0.);
			auto && tmp_data_cache = make_solver_vector_default_wrap(*this, this->data_cache);
			tmp_data_cache.resize();
			tmp_data_cache.clear();

			this->assembleToGlobalArray(dof_id, x, data_cache, 1.);

			auto && tmp_array = make_solver_vector_default_wrap(*this, array);
			tmp_array.resize();
			tmp_array.clear();

			A.matVecMul(tmp_data_cache, tmp_array, scale_factor, 1.);
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::assembleMatMulVectToArray(const ID & dof_id,
			const ID & A_id,
			const Array<Real> & x,
			Array<Real> & array,
			Real scale_factor) {
			Array<Real> tmp_array(0, 1, id + ":tmp_array");
			assembleMatMulVectToGlobalArray(dof_id, A_id, x, tmp_array, scale_factor);
			getArrayPerDOFs(dof_id, tmp_array, array);
			}


			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::assembleMatMulVectToResidual(const ID & dof_id,
			const ID & A_id,
			const Array<Real> & x,
			Real scale_factor) {
			Array<Real> tmp_residual(0, 1, id + ":tmp_residual");
			assembleMatMulVectToGlobalArray(dof_id, A_id, x, tmp_residual, scale_factor);
			this->getResidualArray() += tmp_residual;
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::assembleLumpedMatMulVectToResidual(
			const ID & dof_id, const ID & A_id, const Array<Real> & x,
			Real scale_factor) {
			const Array<Real> & A = this->getLumpedMatrix(A_id);

			this->data_cache.resize(this->local_system_size);
			this->data_cache.clear();
			this->assembleToGlobalArray(dof_id, x, data_cache, scale_factor);

			for (auto && data : zip(make_view(A), make_view(data_cache),
			make_view(this->getResidualArray()))) {
			const auto & A = std::get<0>(data);
			const auto & x = std::get<1>(data);
			auto & r = std::get<2>(data);
			r += A * x;
			}
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::assembleElementalMatricesToMatrix(
			const ID & matrix_id, const ID & dof_id, const Array<Real> & elementary_mat,
			const ElementType & type, const GhostType & ghost_type,
			const MatrixType & elemental_matrix_type,
			const Array<UInt> & filter_elements) {
			AKANTU_DEBUG_IN();

			auto & dof_data = this->getDOFData(dof_id);

			AKANTU_DEBUG_ASSERT(dof_data.support_type == _dst_nodal,
			"This function applies only on Nodal dofs");

			this->addToProfile(matrix_id, dof_id, type, ghost_type);

			const auto & equation_number = this->getLocalEquationsNumbers(dof_id);
			auto & A = this->getMatrix(matrix_id);

			UInt nb_element;
			UInt * filter_it = nullptr;
			if (filter_elements != empty_filter) {
			nb_element = filter_elements.size();
			filter_it = filter_elements.storage();
			} else {
			if (dof_data.group_support != "__mesh__") {
			const auto & group_elements =
			this->mesh->getElementGroup(dof_data.group_support)
			.getElements(type, ghost_type);
			nb_element = group_elements.size();
			filter_it = group_elements.storage();
			} else {
			nb_element = this->mesh->getNbElement(type, ghost_type);
			}
			}

			AKANTU_DEBUG_ASSERT(elementary_mat.size() == nb_element,
			"The vector elementary_mat("
			<< elementary_mat.getID()
			<< ") has not the good size.");

			UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);

			UInt nb_degree_of_freedom = dof_data.dof->getNbComponent();

			const Array<UInt> & connectivity =
			this->mesh->getConnectivity(type, ghost_type);
			auto conn_begin = connectivity.begin(nb_nodes_per_element);
			auto conn_it = conn_begin;
			auto size_mat = nb_nodes_per_element * nb_degree_of_freedom;

			Vector<Int> element_eq_nb(nb_degree_of_freedom * nb_nodes_per_element);
			auto el_mat_it = elementary_mat.begin(size_mat, size_mat);

			for (UInt e = 0; e < nb_element; ++e, ++el_mat_it) {
			if (filter_it)
			conn_it = conn_begin + *filter_it;

			this->extractElementEquationNumber(equation_number, *conn_it,
			nb_degree_of_freedom, element_eq_nb);
			std::transform(element_eq_nb.begin(), element_eq_nb.end(),
			element_eq_nb.begin(), [&](auto && local) {
			return this->localToGlobalEquationNumber(local);
			});

			if (filter_it)
			++filter_it;
			else
			++conn_it;

			if (A.getMatrixType() == _symmetric)
			if (elemental_matrix_type == _symmetric)
			this->addSymmetricElementalMatrixToSymmetric(A, *el_mat_it,
			element_eq_nb, A.size());
			else
			this->addUnsymmetricElementalMatrixToSymmetric(A, *el_mat_it,
			element_eq_nb, A.size());
			else
			this->addElementalMatrixToUnsymmetric(A, *el_mat_it, element_eq_nb,
			A.size());
			}

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::assemblePreassembledMatrix(
			const ID & dof_id_m, const ID & dof_id_n, const ID & matrix_id,
			const TermsToAssemble & terms) {
			const auto & equation_number_m = this->getLocalEquationsNumbers(dof_id_m);
			const auto & equation_number_n = this->getLocalEquationsNumbers(dof_id_n);
			SparseMatrixAIJ & A = this->getMatrix(matrix_id);

			for (const auto & term : terms) {
			auto gi = this->localToGlobalEquationNumber(equation_number_m(term.i()));
			auto gj = this->localToGlobalEquationNumber(equation_number_n(term.j()));
			A.add(gi, gj, term);
			}
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::addToProfile(const ID & matrix_id, const ID & dof_id,
			const ElementType & type,
			const GhostType & ghost_type) {
			AKANTU_DEBUG_IN();

			const auto & dof_data = this->getDOFData(dof_id);

			if (dof_data.support_type != _dst_nodal)
			return;

			auto mat_dof = std::make_pair(matrix_id, dof_id);
			auto type_pair = std::make_pair(type, ghost_type);

			auto prof_it = this->matrix_profiled_dofs.find(mat_dof);
			if (prof_it != this->matrix_profiled_dofs.end() &&
			std::find(prof_it->second.begin(), prof_it->second.end(), type_pair) !=
			prof_it->second.end())
			return;

			auto nb_degree_of_freedom_per_node = dof_data.dof->getNbComponent();

			const auto & equation_number = this->getLocalEquationsNumbers(dof_id);

			auto & A = this->getMatrix(matrix_id);

			auto size = A.size();

			auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);

			const auto & connectivity = this->mesh->getConnectivity(type, ghost_type);
			auto cbegin = connectivity.begin(nb_nodes_per_element);
			auto cit = cbegin;

			auto nb_elements = connectivity.size();
			UInt * ge_it = nullptr;
			if (dof_data.group_support != "__mesh__") {
			const auto & group_elements =
			this->mesh->getElementGroup(dof_data.group_support)
			.getElements(type, ghost_type);
			ge_it = group_elements.storage();
			nb_elements = group_elements.size();
			}

			UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom_per_node;
			Vector<Int> element_eq_nb(size_mat);

			for (UInt e = 0; e < nb_elements; ++e) {
			if (ge_it)
			cit = cbegin + *ge_it;

			this->extractElementEquationNumber(
			equation_number, *cit, nb_degree_of_freedom_per_node, element_eq_nb);
			std::transform(
			element_eq_nb.storage(), element_eq_nb.storage() + element_eq_nb.size(),
			element_eq_nb.storage(),
			[&](auto & local) { return this->localToGlobalEquationNumber(local); });

			if (ge_it)
			++ge_it;
			else
			++cit;

			for (UInt i = 0; i < size_mat; ++i) {
			UInt c_irn = element_eq_nb(i);
			if (c_irn < size) {
			for (UInt j = 0; j < size_mat; ++j) {
			UInt c_jcn = element_eq_nb(j);
			if (c_jcn < size) {
			A.add(c_irn, c_jcn);
			}
			}
			}
			}
			}

			this->matrix_profiled_dofs[mat_dof].push_back(type_pair);

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::applyBoundary(const ID & matrix_id) {
			SparseMatrixAIJ & J = this->getMatrix(matrix_id);

			if (this->jacobian_release == J.getValueRelease()) {
			auto are_equal =
			std::equal(global_blocked_dofs.begin(), global_blocked_dofs.end(),
			previous_global_blocked_dofs.begin());

			if (not are_equal)
			J.applyBoundary();

			previous_global_blocked_dofs.copy(global_blocked_dofs);
			} else {
			J.applyBoundary();
			}

			this->jacobian_release = J.getValueRelease();
			}

			/* -------------------------------------------------------------------------- */
			// const Array<Real> & DOFManagerDefault::getGlobalResidual() {
			// if (this->synchronizer) {
			// if (not this->global_residual) {
			// this->global_residual =
			// std::make_unique<Array<Real>>(0, 1, "global_residual");
			// }

			// if (this->synchronizer->getCommunicator().whoAmI() == 0) {
			// this->global_residual->resize(this->system_size);
			// this->synchronizer->gather(this->residual, *this->global_residual);
			// } else {
			// this->synchronizer->gather(this->residual);
			// }

			// return *this->global_residual;
			// } else {
			// return this->residual;
			// }
			// }

			/* -------------------------------------------------------------------------- */
			Array<Real> & DOFManagerDefault::getSolutionArray() {
			return dynamic_cast<SolverVectorDefault *>(this->solution.get())->getVector();
			}

			/* -------------------------------------------------------------------------- */
			const Array<Real> & DOFManagerDefault::getResidualArray() const {
			return dynamic_cast<SolverVectorDefault *>(this->residual.get())->getVector();
			}

			/* -------------------------------------------------------------------------- */
			Array<Real> & DOFManagerDefault::getResidualArray() {
			return dynamic_cast<SolverVectorDefault *>(this->residual.get())->getVector();
			}

			/* -------------------------------------------------------------------------- */
			// void DOFManagerDefault::setGlobalSolution(const Array<Real> & solution) {
			// if (this->synchronizer) {
			// if (this->synchronizer->getCommunicator().whoAmI() == 0) {
			// this->synchronizer->scatter(this->global_solution, solution);
			// } else {
			// this->synchronizer->scatter(this->global_solution);
			// }
			// } else {
			// AKANTU_DEBUG_ASSERT(solution.size() == this->global_solution.size(),
			// "Sequential call to this function needs the solution
			// " "to be the same size as the global_solution");
			// this->global_solution.copy(solution);
			// }
			// }

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::onNodesAdded(const Array<UInt> & nodes_list,
			const NewNodesEvent & event) {
			DOFManager::onNodesAdded(nodes_list, event);

			if (this->synchronizer)
			this->synchronizer->onNodesAdded(nodes_list);
			}

			/* -------------------------------------------------------------------------- */
			std::pair<UInt, UInt>
			DOFManagerDefault::updateNodalDOFs(const ID & dof_id,
			const Array<UInt> & nodes_list) {
			UInt nb_new_local_dofs, nb_new_pure_local;
			std::tie(nb_new_local_dofs, nb_new_pure_local) =
			DOFManager::updateNodalDOFs(dof_id, nodes_list);

			auto & dof_data = this->getDOFDataTyped<DOFDataDefault>(dof_id);
			updateDOFsData(dof_data, nb_new_local_dofs, nb_new_pure_local,
			nodes_list.size(),
			[&nodes_list](UInt pos) -> UInt { return nodes_list[pos]; });

			return std::make_pair(nb_new_local_dofs, nb_new_pure_local);
			}

			/* -------------------------------------------------------------------------- */
			/* -------------------------------------------------------------------------- */
			class GlobalDOFInfoDataAccessor : public DataAccessor<UInt> {
			public:
			using size_type =
			typename std::unordered_map<UInt, std::vector<UInt>>::size_type;

			GlobalDOFInfoDataAccessor(DOFManagerDefault::DOFDataDefault & dof_data,
			DOFManagerDefault & dof_manager)
			: dof_data(dof_data), dof_manager(dof_manager) {
			for (auto && pair :
			zip(dof_data.local_equation_number, dof_data.associated_nodes)) {
			UInt node, dof;
			std::tie(dof, node) = pair;

			dofs_per_node[node].push_back(dof);
			}
			}

			UInt getNbData(const Array<UInt> & nodes,
			const SynchronizationTag & tag) const override {
			if (tag == SynchronizationTag::_ask_nodes or
			tag == SynchronizationTag::_giu_global_conn) {
			return nodes.size() * dof_data.dof->getNbComponent() * sizeof(UInt);
			}

			return 0;
			}

			void packData(CommunicationBuffer & buffer, const Array<UInt> & nodes,
			const SynchronizationTag & tag) const override {
			if (tag == SynchronizationTag::_ask_nodes or
			tag == SynchronizationTag::_giu_global_conn) {
			for (auto & node : nodes) {
			auto & dofs = dofs_per_node.at(node);
			for (auto & dof : dofs) {
			buffer << dof_manager.global_equation_number(dof);
			}
			}
			}
			}

			void unpackData(CommunicationBuffer & buffer, const Array<UInt> & nodes,
			const SynchronizationTag & tag) override {
			if (tag == SynchronizationTag::_ask_nodes or
			tag == SynchronizationTag::_giu_global_conn) {
			for (auto & node : nodes) {
			auto & dofs = dofs_per_node[node];
			for (auto dof : dofs) {
			UInt global_dof;
			buffer >> global_dof;
			AKANTU_DEBUG_ASSERT(
			(dof_manager.global_equation_number(dof) == UInt(-1) or
			dof_manager.global_equation_number(dof) == global_dof),
			"This dof already had a global_dof_id which is different from "
			"the received one. "
			<< dof_manager.global_equation_number(dof)
			<< " != " << global_dof);
			dof_manager.global_equation_number(dof) = global_dof;
			dof_manager.global_to_local_mapping[global_dof] = dof;
			}
			}
			}
			}

			protected:
			std::unordered_map<UInt, std::vector<UInt>> dofs_per_node;
			DOFManagerDefault::DOFDataDefault & dof_data;
			DOFManagerDefault & dof_manager;
			};

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::resizeGlobalArrays() {
			// resize all relevant arrays
			this->residual->resize();
			this->solution->resize();

			this->dofs_flag.resize(this->local_system_size, NodeFlag::_normal);
			this->global_blocked_dofs.resize(this->local_system_size, true);
			this->previous_global_blocked_dofs.resize(this->local_system_size, true);
			this->global_equation_number.resize(this->local_system_size, -1);

			for (auto & lumped_matrix : lumped_matrices)
			lumped_matrix.second->resize(this->local_system_size);

			matrix_profiled_dofs.clear();
			for (auto & matrix : matrices) {
			matrix.second->clearProfile();
			}
			}

			/* -------------------------------------------------------------------------- */
			auto DOFManagerDefault::computeFirstDOFIDs(UInt nb_new_local_dofs,
			UInt nb_new_pure_local) {
			auto prank = this->communicator.whoAmI();
			auto psize = this->communicator.getNbProc();

			// determine the first local/global dof id to use
			Array<UInt> nb_dofs_per_proc(psize);
			nb_dofs_per_proc(prank) = nb_new_pure_local;
			this->communicator.allGather(nb_dofs_per_proc);

			this->first_global_dof_id += std::accumulate(
			nb_dofs_per_proc.begin(), nb_dofs_per_proc.begin() + prank, 0);

			auto first_global_dof_id = this->first_global_dof_id;
			auto first_local_dof_id = this->local_system_size - nb_new_local_dofs;

			this->first_global_dof_id += std::accumulate(nb_dofs_per_proc.begin() + prank,
			nb_dofs_per_proc.end(), 0);

			return std::make_pair(first_local_dof_id, first_global_dof_id);
			}
			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::updateDOFsData(
			DOFDataDefault & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local,
			UInt nb_node, const std::function<UInt(UInt)> & getNode) {
			auto nb_local_dofs_added = nb_node * dof_data.dof->getNbComponent();

			resizeGlobalArrays();

			auto first_dof_pos = dof_data.local_equation_number.size();
			dof_data.local_equation_number.reserve(dof_data.local_equation_number.size() +
			nb_local_dofs_added);
			dof_data.associated_nodes.reserve(dof_data.associated_nodes.size() +
			nb_local_dofs_added);

			std::unordered_map<std::pair<UInt, UInt>, UInt> masters_dofs;

			// update per dof info
			UInt local_eq_num, first_global_dof_id;
			std::tie(local_eq_num, first_global_dof_id) =
			computeFirstDOFIDs(nb_new_local_dofs, nb_new_pure_local);
			for (auto d : arange(nb_local_dofs_added)) {
			auto node = getNode(d / dof_data.dof->getNbComponent());
			auto dof_flag = this->mesh->getNodeFlag(node);

			dof_data.associated_nodes.push_back(node);
			auto is_local_dof = this->mesh->isLocalOrMasterNode(node);
			auto is_periodic_slave = this->mesh->isPeriodicSlave(node);
			auto is_periodic_master = this->mesh->isPeriodicMaster(node);

			if (is_periodic_slave) {
			dof_data.local_equation_number.push_back(UInt(-1));
			continue;
			}

			// update equation numbers
			this->dofs_flag(local_eq_num) = dof_flag;
			dof_data.local_equation_number.push_back(local_eq_num);

			if (is_local_dof) {
			this->global_equation_number(local_eq_num) = first_global_dof_id;
			this->global_to_local_mapping[first_global_dof_id] = local_eq_num;
			++first_global_dof_id;
			} else {
			this->global_equation_number(local_eq_num) = UInt(-1);
			}

			if (is_periodic_master) {
			auto node = getNode(d / dof_data.dof->getNbComponent());
			auto dof = d % dof_data.dof->getNbComponent();
			masters_dofs.insert(
			std::make_pair(std::make_pair(node, dof), local_eq_num));
			}

			++local_eq_num;
			}

			// correct periodic slave equation numbers
			if (this->mesh->isPeriodic()) {
			auto assoc_begin = dof_data.associated_nodes.begin();
			for (auto d : arange(nb_local_dofs_added)) {
			auto node = dof_data.associated_nodes(first_dof_pos + d);
			if (not this->mesh->isPeriodicSlave(node))
			continue;

			auto master_node = this->mesh->getPeriodicMaster(node);
			auto dof = d % dof_data.dof->getNbComponent();
			dof_data.local_equation_number(first_dof_pos + d) =
			masters_dofs[std::make_pair(master_node, dof)];
			}
			}

			// synchronize the global numbering for slaves
			if (this->synchronizer) {
			GlobalDOFInfoDataAccessor data_accessor(dof_data, *this);

			if (this->mesh->isPeriodic()) {
			mesh->getPeriodicNodeSynchronizer().synchronizeOnce(
			data_accessor, SynchronizationTag::_giu_global_conn);
			}

			auto & node_synchronizer = this->mesh->getNodeSynchronizer();
			node_synchronizer.synchronizeOnce(data_accessor,
			SynchronizationTag::_ask_nodes);
			}
			}

			/* -------------------------------------------------------------------------- */
			void DOFManagerDefault::updateDOFsData(DOFDataDefault & dof_data,
			UInt nb_new_local_dofs,
			UInt nb_new_pure_local) {
			resizeGlobalArrays();

			dof_data.local_equation_number.reserve(dof_data.local_equation_number.size() +
			nb_new_local_dofs);

			UInt first_local_dof_id, first_global_dof_id;
			std::tie(first_local_dof_id, first_global_dof_id) =
			computeFirstDOFIDs(nb_new_local_dofs, nb_new_pure_local);

			// update per dof info
			for (auto _ [[gnu::unused]] : arange(nb_new_local_dofs)) {
			// update equation numbers
			this->dofs_flag(first_local_dof_id) = NodeFlag::_normal;
			;
			dof_data.local_equation_number.push_back(first_local_dof_id);

			this->global_equation_number(first_local_dof_id) = first_global_dof_id;
			this->global_to_local_mapping[first_global_dof_id] = first_local_dof_id;

			++first_global_dof_id;
			++first_local_dof_id;
			}
			}

			/* -------------------------------------------------------------------------- */
			// register in factory
			static bool default_dof_manager_is_registered [[gnu::unused]] =
			DefaultDOFManagerFactory::getInstance().registerAllocator(
			"default",
			[](const ID & id,
			const MemoryID & mem_id) -> std::unique_ptr<DOFManager> {
			return std::make_unique<DOFManagerDefault>(id, mem_id);
			});

			static bool dof_manager_is_registered [[gnu::unused]] =
			DOFManagerFactory::getInstance().registerAllocator(
			"default",
			[](Mesh & mesh, const ID & id,
			const MemoryID & mem_id) -> std::unique_ptr<DOFManager> {
			return std::make_unique<DOFManagerDefault>(mesh, id, mem_id);
			});
			} // namespace akantu

rAKA/src/model/common/dof_manager/dof_manager_default.cc05a352d13d99features/34-poisson…

dof_manager_default.ccView Options

rAKA/src/model/common/dof_manager/dof_manager_default.cc
05a352d13d99features/34-poisson…

dof_manager_default.cc
View Options