shape_matrix.hh
3a940948c399master

shape_matrix.hh
View Options

			/**
			* @file shape_matrix.hh
			*
			* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
			*
			* @date Fri Jul 11 15:47:44 2014
			*
			* @brief Data structure for the storage of shape function values on
			* atomic/point positions
			*
			* @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.
			*
			*/

			#ifndef __LIBMULTISCALE_SHAPE_MATRIX_HH__
			#define __LIBMULTISCALE_SHAPE_MATRIX_HH__
			/* -------------------------------------------------------------------------- */
			#include "container_array.hh"
			/* ------------------------------------------------------------------------ */
			#include <Eigen/Sparse>
			/* ------------------------------------------------------------------------ */

			__BEGIN_LIBMULTISCALE__

			/* ------------------------------------------------------------------------ */

			template <typename T> struct SwapableObject {

			SwapableObject() {
			current = &obj1;
			_next = &obj2;
			}

			void swap() {
			auto *temp = current;
			current = _next;
			_next = temp;
			}
			T &operator() { return current; };
			T *operator->() { return current; };
			template <typename T1> decltype(auto) operator[](T1 &&idx) {
			return (*current)[idx];
			}

			T &next() { return *_next; }

			T obj1;
			T obj2;
			T *current;
			T *_next;
			};

			/* ------------------------------------------------------------------------ */

			/**
			Class ShapeMatrix

			This class is actually a sparsematrix to store the
			shape functions of several mesh nodes evaluated on
			specific locations.

			The structure follows:

			$m_{Ii} = \varphi_I(X_i)$

			where $I$ is a node index, and $i$ an atom index.

			*/

			class ShapeMatrix
			: public Eigen::SparseMatrix<libmultiscale::Real, Eigen::RowMajor> {
			public:
			using Eigen::SparseMatrix<libmultiscale::Real, Eigen::RowMajor>::SparseMatrix;
			~ShapeMatrix();

			inline Real interpolate(ContainerArray<Real> &atom_data,
			ContainerArray<Real> &node_data);
			inline void removeAtom(UInt at, UInt elem);
			inline void changeAtomIndirection(UInt el, UInt old_ind_atom,
			UInt new_ind_atom);
			inline void swapAtomIndirections();

			//! return true if the matrix is ready to be used
			// bool isFilled();

			private:
			// inline void setSub(UInt i, UInt m, UInt n);
			// inline void allocate();
			// inline void print();
			// inline void fill(UInt atom, UInt el, std::vector<UInt> &node_global,
			// std::vector<UInt> &nodes, std::vector<Real> &values);
			// template <UInt Dim, typename _Vec_>
			// inline Vector<Dim> interpolate(UInt at, _Vec_ &field, UInt elem);

			inline void getRelatedAtomsToElem(UInt elem, std::vector<UInt> &atoms);
			inline void getRelatedNodesToElem(UInt elem, std::vector<UInt> &nodes);
			inline void averageOnElements(ContainerArray<Real> &DataA,
			ContainerArray<Real> &DataM);

			inline void
			alterMatrixIndexesForPBC(std::vector<std::pair<UInt, UInt>> &pbc_pairs);

			public:
			/* ------------------------------------------------------------------------ */
			/* Methods (specific to a coupling type) */
			/* ------------------------------------------------------------------------ */

			inline void buildNodeShape(ContainerArray<Real> &node_shape);
			inline void buildThermalStuff(std::vector<Real> &nodal_out,
			std::vector<Real> &atom_shape);
			inline void buildContinuumConstraint(ContainerArray<Real> &A,
			ContainerArray<Real> &node_shape,
			ContainerArray<Real> &node_weight);
			template <UInt Dim, typename Vec>
			void buildRHS(ContainerArray<Real> &rhs, Vec &v);
			template <UInt Dim, typename _Vec_>
			inline void applyContinuumConstraint(_Vec_ &v, ContainerArray<Real> &multL,
			ContainerArray<Real> &lambdas);
			template <UInt Dim, typename _Vec_>
			inline void applyChangingThermal(_Vec_ &vDisp, _Vec_ &vVel, _Vec_ &vAcc,
			_Vec_ &masses, std::vector<Real> &multL,
			Real deltat);
			inline void buildLeastSquareMatrix(Eigen::MatrixXf &A);
			inline void buildLeastSquareRHS(ContainerArray<Real> &rhs,
			ContainerArray<Real> &v);
			inline void buildFormeFaibleCMatrix(Eigen::MatrixXf &C,
			std::vector<Real> &lambdas);
			inline void correctFormeFaible(Real &v, Array &rhs, std::vector<Real> &lamb,
			UInt at, UInt elem);

			inline auto indirectionAtomTable();
			inline void extractShapeMatrix(Eigen::MatrixXf &A);

			private:
			//! filled indicator
			// bool is_filled = false;
			//! number of nodes in the mesh
			// UInt nb_nodes;
			//! number of points
			// UInt nb_points;
			//! number of stored elements
			// UInt nb_elements;
			//! size of submatrices (m1,n1,m2,n2,...)
			// std::vector<std::pair<UInt, UInt>> matrix_sizes;
			//! array of submatrices
			// std::vector<Eigen::MatrixXd> SubMatrix;
			//! map element index to contained atom indexes
			// SwapableObject<std::vector<std::vector<UInt>>> elementToAtoms;
			//! map element index to contained node indexes
			// std::vector<std::vector<UInt>> elementToNodes;
			//! index map to node in the global index (big vector)
			// std::vector<std::vector<UInt>> elementToNodesGlobalIndex;
			};

			/* -------------------------------------------------------------------------- */
			// inline ShapeMatrix::ShapeMatrix(UInt nb_elements, UInt nb_nodes,
			// UInt nb_points) {
			// this->nb_elements = nb_elements;
			// this->nb_points = nb_points;
			// this->nb_nodes = nb_nodes;
			// DUMP("creating shapematrix for " << nb_elements << " elements",
			// DBG_WARNING); SubMatrix.resize(nb_elements);
			// matrix_sizes.resize(nb_elements);
			// elementToAtoms.obj1.resize(nb_elements);
			// elementToAtoms.obj2.resize(nb_elements);
			// elementToNodes.resize(nb_elements);
			// elementToNodesGlobalIndex.resize(nb_elements);
			// }
			/* -------------------------------------------------------------------------- */

			// inline ShapeMatrix::ShapeMatrix(UInt nb_nodes, UInt nb_points) {
			// this->nb_points = nb_points;
			// this->nb_nodes = nb_nodes;
			// DUMP("creating shapematrix for " << nb_elements << " elements",
			// DBG_WARNING); SubMatrix.resize(nb_elements);
			// matrix_sizes.resize(nb_elements);
			// elementToAtoms.obj1.resize(nb_elements);
			// elementToAtoms.obj2.resize(nb_elements);
			// elementToNodes.resize(nb_elements);
			// elementToNodesGlobalIndex.resize(nb_elements);
			// }
			// /* --------------------------------------------------------------------------
			// */

			// inline ShapeMatrix::ShapeMatrix(UInt nb_elements) {
			// LM_TOIMPLEMENT;
			// this->nb_elements = nb_elements;
			// DUMP("creating shapematrix for " << nb_elements << " elements",
			// DBG_WARNING); SubMatrix.resize(nb_elements);
			// matrix_sizes.resize(nb_elements);
			// elementToAtoms.obj1.resize(nb_elements);
			// elementToAtoms.obj2.resize(nb_elements);
			// elementToNodes.resize(nb_elements);
			// elementToNodesGlobalIndex.resize(nb_elements);
			// }
			/* -------------------------------------------------------------------------- */
			inline ShapeMatrix::~ShapeMatrix() {}
			/* -------------------------------------------------------------------------- */
			// inline bool ShapeMatrix::isFilled() { return this->is_filled; }
			/* -------------------------------------------------------------------------- */

			// inline void ShapeMatrix::setSub(UInt i, UInt m, UInt n) {
			// LM_ASSERT(m > 0 && n > 0,
			// "can not set a smatrix bloc size of zero : m=" << m << " n=" <<
			// n);
			// matrix_sizes[i].first = m;
			// matrix_sizes[i].second = n;
			// }
			/* -------------------------------------------------------------------------- */
			// inline void ShapeMatrix::allocate() {
			// // first pass to compute the global ShapeMatrix coefficients

			// UInt m, n;

			// UInt total_atoms = 0;
			// UInt total_nodes = 0;

			// for (UInt i = 0; i < nb_elements; ++i) {
			// total_atoms += matrix_sizes[i].first;
			// total_nodes += matrix_sizes[i].second;
			// }

			// UInt index_atoms = 0;
			// UInt index_nodes = 0;

			// // second pass to allocate the data structures of matrices
			// for (UInt i = 0; i < nb_elements; ++i) {
			// m = matrix_sizes[i].first;
			// n = matrix_sizes[i].second;

			// SubMatrix[i].resize(m + 1, n);
			// SubMatrix[i].setZero();

			// // increment of indexes
			// index_atoms += m;
			// index_nodes += n;
			// }
			// }
			/* -------------------------------------------------------------------------- */
			// inline void ShapeMatrix::print() {
			// std::cout << *this << std::endl;
			// for (UInt k = 0; k < nb_elements; ++k) {
			// auto &mat = SubMatrix[k];
			// if (mat.rows() == 1)
			// return;

			// DUMP("affiche les infos de " << k, DBG_DETAIL);

			// for (UInt i = 0; i < mat.rows() - 1; ++i)
			// std::cout << elementToAtoms[k][i] << "\t";

			// std::cout << std::endl;

			// std::cout << SubMatrix[k] << std::endl;
			// }
			//}
			/* -------------------------------------------------------------------------- */
			inline void ShapeMatrix::changeAtomIndirection(UInt el, UInt old_ind_atom,
			UInt new_ind_atom) {
			LM_TOIMPLEMENT;
			// needs to swap column old_ind_atom with new_ind_atom
			auto c1 = this->col(old_ind_atom);
			auto c2 = this->col(new_ind_atom);
			H.col(i) = Hj;
			H.col(j) = Hi;

			// auto &mat = SubMatrix[el];
			// auto &indA = elementToAtoms[el];
			// auto &indAN = elementToAtoms.next()[el];
			// indAN.resize(indA.size());

			// // should be done in swapatomindirection TODO !!
			// for (UInt i = 0; i < mat.rows() - 1; ++i)
			// if (indA[i] == UINT_MAX)
			// indAN[i] = UINT_MAX;

			// UInt i = 0;
			// for (i = 0; i < mat.rows() - 1; ++i) {
			// if (indA[i] == old_ind_atom) {
			// indAN[i] = new_ind_atom;
			// break;
			// }
			// }
			// LM_ASSERT(i != mat.rows() - 1,
			// "atom not found this should not append (index was " <<
			// old_ind_atom
			// << ")");
			}
			/* --------------------------------------------------------------------------
			*/
			// inline void ShapeMatrix::fill(UInt atom, UInt el,
			// std::vector<UInt> &nodes_global,
			// std::vector<UInt> &nodes,
			// std::vector<Real> &values) {

			// UInt atom_index = elementToAtoms[el].size();
			// elementToAtoms[el].push_back(atom);

			// if (elementToNodes[el].size() == 0) {
			// elementToNodes[el] = nodes;
			// elementToNodesGlobalIndex[el] = nodes_global;
			// }

			// for (UInt nd = 0; nd < nodes.size(); ++nd) {

			// Real &value = values[nd];
			// UInt &node_global [[gnu::unused]] = nodes_global[nd];
			// UInt &node [[gnu::unused]] = nodes[nd];
			// LM_ASSERT(!std::isnan(value) && !std::isinf(value),
			// "the passed value for atom " << atom << " node_global "
			// << node_global << " node " << node
			// << "el " << el);

			// DUMP("setting shapes for atom " << atom << " and element " << el << "
			// node "
			// << node << " node_global " << node_global
			// << " value " << value,
			// DBG_ALL);
			// auto &m = SubMatrix[el];

			// m(atom_index, nd) = value;
			// // make the sum shapes for all atoms
			// m(m.rows() - 1, nd) += value;
			// }
			// }
			/* --------------------------------------------------------------------------
			*/
			inline void ShapeMatrix::removeAtom(UInt at, UInt elem) {

			LM_TOIMPLEMENT;
			// Do not really remove the atom.
			// Only substract the contribution to the cumulative part in submatrix.
			auto &mat = SubMatrix[elem];
			auto &indA = elementToAtoms[elem];
			UInt i = 0;
			for (i = 0; i < mat.rows() - 1; ++i)
			if (indA[i] == at)
			break;

			LM_ASSERT(i != mat.rows() - 1, "could not find the atom to remove !!");

			for (UInt j = 0; j < mat.cols(); ++j) {
			mat(mat.rows() - 1, j) -= mat(i, j);
			}
			indA[i] = UINT_MAX;
			}
			/* --------------------------------------------------------------------------
			*/
			inline void ShapeMatrix::buildNodeShape(ContainerArray<Real> &node_shape) {
			for (UInt k = 0; k < nb_elements; ++k) {
			auto &indN = elementToNodes[k];
			auto &mat = SubMatrix[k];
			for (UInt j = 0; j < mat.cols(); ++j) {
			DUMP("Summing node_shape[" << indN[j] << " (el= " << k << ", j = " << j
			<< ")] with " << mat(mat.rows() - 1, j),
			DBG_ALL);
			LM_ASSERT(mat(mat.rows() - 1, j) >= 0,
			"wrong shape fonction in m("
			<< mat.rows() - 1 << ", " << j
			<< ") = " << &mat(mat.rows() - 1, j) << " = "
			<< mat(mat.rows() - 1, j));
			node_shape[indN[j]] += mat(mat.rows() - 1, j);
			}
			}
			}
			/* --------------------------------------------------------------------------
			*/
			inline void
			ShapeMatrix::buildThermalStuff([[gnu::unused]] std::vector<Real> &nodal_out,
			[[gnu::unused]] std::vector<Real> &atom_shape) {
			LM_TOIMPLEMENT;
			// for (UInt elem = 0; elem < nb_elements; ++elem) {
			// UInt *indN = elementToNodes[elem];
			// auto &indA = elementToAtoms[elem];
			// auto &mat = SubMatrix[elem];
			// for (UInt i = 0; i < mat.rows() - 1; ++i) {
			// for (UInt j = 0; j < mat.cols(); ++j) {
			// nodal_out[indN[j]] += mat(i, j) * atom_shape[indA[i]];
			// }
			// }
			// }
			}
			/* --------------------------------------------------------------------------
			*/
			inline void
			ShapeMatrix::buildContinuumConstraint(ContainerArray<Real> &A,
			ContainerArray<Real> &node_shape,
			ContainerArray<Real> &node_weight) {
			for (UInt k = 0; k < nb_elements; ++k) {
			auto &indA = elementToAtoms[k];
			auto &indN = elementToNodes[k];
			auto &mat = SubMatrix[k];
			for (UInt i = 0; i < mat.rows() - 1; ++i)
			if (indA[i] != UINT_MAX)
			for (UInt j = 0; j < mat.cols(); ++j) {
			DUMP("Building A["
			<< indA[i] << "]"
			<< " += " << mat(i, j) << "*"
			<< " " << node_shape[indN[j]] << " / "
			<< node_weight[indN[j]] << std::endl
			<< " = "
			<< mat(i, j) * node_shape[indN[j]] / node_weight[indN[j]],
			DBG_ALL);
			A[indA[i]] += mat(i, j) * node_shape[indN[j]] / node_weight[indN[j]];
			LM_ASSERT(!std::isinf(A[indA[i]]) && !std::isnan(A[indA[i]]),
			"node_weight[" << indN[j] << "] = " << node_weight[indN[j]]
			<< "node_shape[" << indN[j]
			<< "] = " << node_shape[indN[j]]
			<< "fail on atom " << i);
			}
			}
			}
			/* --------------------------------------------------------------------------
			*/
			// template <UInt Dim, typename _Vec_>
			// void ShapeMatrix::buildRHS(ContainerArray<Real> &rhs, _Vec_ &v) {

			// for (UInt k = 0; k < nb_elements; ++k) {
			// auto &mat = SubMatrix[k];
			// auto &indGN = elementToNodesGlobalIndex[k];
			// auto &indA = elementToAtoms[k];
			// // loop over all the atoms of the bloc
			// for (UInt i = 0; i < mat.rows() - 1; ++i)
			// if (indA[i] != UINT_MAX) {
			// VectorView<Dim> _rhs(&rhs(indA[i], 0));
			// // summing contributions ovre the nodes
			// for (UInt j = 0; j < mat.cols(); ++j) {
			// VectorView<Dim> _v(&v(indGN[j], 0));
			// DUMP("Before rhs[" << indA[i] << "]=" << _rhs, DBG_ALL);
			// _rhs += _v * mat(i, j);
			// DUMP("Building " << _rhs << "=rhs[" << indA[i] << "] with "
			// << mat(i, j) << "" << _v << " = " << _v
			// mat(i, j)
			// << " (index = " << indGN[j] << ")",
			// DBG_ALL);
			// }
			// }
			// }
			// }
			/* --------------------------------------------------------------------------
			*/
			template <UInt Dim, typename _Vec_>
			void ShapeMatrix::applyContinuumConstraint(_Vec_ &velocity,
			ContainerArray<Real> &multL,
			ContainerArray<Real> &lambdas) {
			// loop over the blocs
			for (UInt k = 0; k < nb_elements; ++k) {
			auto &mat = SubMatrix[k];
			// get indirection to global node indices
			// get direct access to array
			auto &indGN = elementToNodesGlobalIndex[k];
			auto &indN = elementToNodes[k];
			auto &indA = elementToAtoms[k];
			// loop over the atoms of the bloc
			for (UInt i = 0; i < mat.rows() - 1; ++i)
			if (indA[i] != UINT_MAX) {
			VectorView<Dim> _multL(&multL(indA[i], 0));
			// loop over the nodes
			for (UInt j = 0; j < mat.cols(); ++j) {
			VectorView<Dim> _v(&velocity(indGN[j], 0));
			DUMP("adding v[" << indGN[j] << "] with " << -1.0 * mat(i, j) << " * "
			<< _multL << " {rhs[" << indA[i] << "]}/ "
			<< lambdas[indN[j]] << " = "
			<< -1.0 * mat(i, j) * _multL / lambdas[indN[j]],
			DBG_ALL);

			_v += -1.0 * mat(i, j) * _multL / lambdas[indN[j]];
			}
			}
			}
			}
			/* --------------------------------------------------------------------------
			*/
			template <UInt Dim, typename _Vec_>
			void ShapeMatrix::applyChangingThermal(_Vec_ &vDisp, _Vec_ &vVel, _Vec_ &vAcc,
			_Vec_ &masses, std::vector<Real> &multL,
			Real) {
			// loop over the blocs
			for (UInt k = 0; k < nb_elements; ++k) {
			auto &mat = SubMatrix[k];
			// get indirection to global node indices
			// get direct access to array
			auto &indGN = elementToNodesGlobalIndex[k];
			auto &indA = elementToAtoms[k];
			// loop over the atoms of the block
			for (UInt i = 0; i < mat.rows() - 1; ++i)
			if (indA[i] != UINT_MAX)
			// loop over the nodes
			for (UInt j = 0; j < mat.cols(); ++j) {
			vDisp.add(indGN[j] * Dim, -1.0 * mat(i, j) * multL[indA[i] * Dim] /
			masses[indGN[j] * Dim]);
			vDisp.add(indGN[j] * Dim + 1, -1.0 * mat(i, j) *
			multL[indA[i] * Dim + 1] /
			masses[indGN[j] * Dim + 1]);
			vDisp.add(indGN[j] * Dim + 2, -1.0 * mat(i, j) *
			multL[indA[i] * Dim + 2] /
			masses[indGN[j] * Dim + 2]);

			DUMP("Dispx( " << indGN[j] << "," << vDisp[indGN[j] * Dim]
			<< ")+=" << mat(i, j) << "" << multL[indA[i] Dim]
			<< "/" << masses[indGN[j] * Dim],
			DBG_ALL);
			DUMP("Dispy( " << indGN[j] << "," << vDisp[indGN[j] * Dim + 1]
			<< ")+=" << mat(i, j) << "*"
			<< multL[indA[i] * Dim + 1] << "/"
			<< masses[indGN[j] * Dim + 1],
			DBG_ALL);
			DUMP("Dispz( " << indGN[j] << "," << vDisp[indGN[j] * Dim + 2]
			<< ")+=" << mat(i, j) << "*"
			<< multL[indA[i] * Dim + 2] << "/"
			<< masses[indGN[j] * Dim + 2],
			DBG_ALL);

			vVel(indGN[j] * Dim) = 0;
			vVel(indGN[j] * Dim + 1) = 0;
			vVel(indGN[j] * Dim + 2) = 0;

			vAcc(indGN[j] * Dim) = 0;
			vAcc(indGN[j] * Dim + 1) = 0;
			vAcc(indGN[j] * Dim + 2) = 0;
			}
			}
			}
			/* ------------------------------------------------------------------------ */
			// template <UInt Dim, typename _Vec_>
			// Vector<Dim> ShapeMatrix::interpolate(UInt at, _Vec_ &field, UInt elem) {
			// Vector<Dim> res;
			// res.setZero();
			// LM_ASSERT(elem < nb_elements, "overflow detected");
			// auto &indA = elementToAtoms[elem];
			// auto &indGN = elementToNodesGlobalIndex[elem];
			// auto &mat = SubMatrix[elem];

			// DUMP("interpolating for elem " << elem << " atom " << at, DBG_ALL);
			// UInt i;
			// for (i = 0; i < mat.rows() - 1; ++i)
			// if (indA[i] == at)
			// break;
			// for (UInt j = 0; j < mat.cols(); ++j) {
			// for (UInt deg = 0; deg < Dim; ++deg) {
			// res[deg] += mat(i, j) * field(indGN[j], deg);
			// }
			// }
			// return res;
			// }
			/* --------------------------------------------------------------------------
			*/
			inline Real ShapeMatrix::interpolate(ContainerArray<Real> &atom_data,
			ContainerArray<Real> &node_data) {

			UInt stride = atom_data.cols();
			for (UInt elem = 0; elem < nb_elements; ++elem) {
			auto &indN = elementToNodes[elem];
			auto &indA = elementToAtoms[elem];
			auto &mat = SubMatrix[elem];
			for (UInt i = 0; i < mat.rows() - 1; ++i) {
			for (UInt j = 0; j < mat.cols(); ++j) {
			for (UInt k = 0; k < stride; ++k) {
			atom_data(indA[i], k) += mat(i, j) * node_data(indN[j], k);
			}
			}
			}
			}
			return 0;
			}
			/* --------------------------------------------------------------------------
			*/
			inline void ShapeMatrix::getRelatedAtomsToElem(UInt elem,
			std::vector<UInt> &atoms) {

			auto &indA = elementToAtoms[elem];
			auto &mat = SubMatrix[elem];

			atoms.clear();

			// loop over the atoms
			UInt i;
			for (i = 0; i < mat.rows() - 1; ++i)
			if (indA[i] != UInt(-1))
			atoms.push_back(indA[i]);
			}
			/* --------------------------------------------------------------------------
			*/
			inline void ShapeMatrix::getRelatedNodesToElem(UInt elem,
			std::vector<UInt> &nodes) {

			auto &indGN = elementToNodesGlobalIndex[elem];
			// auto &mat = SubMatrix[elem];

			nodes = indGN;
			}
			/* --------------------------------------------------------------------------
			*/
			inline void ShapeMatrix::extractShapeMatrix(Eigen::MatrixXf &A) {
			A.setZero();
			// loop over elements
			for (UInt k = 0; k < nb_elements; ++k) {
			auto &indA = elementToAtoms[k];
			auto &indN = elementToNodes[k];
			auto &mat = SubMatrix[k];
			// loop over nodes
			for (UInt j = 0; j < mat.cols(); ++j) {
			// loop over atoms
			for (UInt i = 0; i < mat.rows() - 1; ++i) {
			A(indA[i], indN[j]) += mat(i, j);
			}
			}
			}
			}
			/* --------------------------------------------------------------------------
			*/

			inline void ShapeMatrix::buildLeastSquareMatrix(Eigen::MatrixXf &A) {
			A.setZero();

			/* on veut construire les A(h,k) */

			/* le produit de deux shape sur la position d'un atome
			n'est nul que si les deux shapes font partie
			d'un element commun : on a une seul boucle pricipale sur
			les elements */

			for (UInt k = 0; k < nb_elements; ++k) {
			auto &indN = elementToNodes[k];
			auto &mat = SubMatrix[k];

			/* ici on met la Real boucle sur les noeuds pour calculer
			la sous matrice de la grande */

			for (UInt j1 = 0; j1 < mat.cols(); ++j1)
			for (UInt j2 = 0; j2 < mat.cols(); ++j2) {
			/* ici on met la boucle de la somme sur les atomes */
			for (UInt i = 0; i < mat.rows() - 1; ++i) {
			A(indN[j1], indN[j2]) += mat(i, j1) * mat(i, j2);
			DUMP("A[" << indN[j1] << "," << indN[j2]
			<< "] = " << A(indN[j1], indN[j2]),
			DBG_ALL);
			}
			}
			}

			std::ofstream f("least-square-constraint.mat");
			f << A;
			}
			/* --------------------------------------------------------------------------
			*/
			inline void ShapeMatrix::buildLeastSquareRHS(ContainerArray<Real> &rhs,
			ContainerArray<Real> &v) {

			UInt stride = rhs.cols();
			DUMP("nb_elements: " << nb_elements, DBG_DETAIL);
			// loop on elements
			for (UInt elem = 0; elem < nb_elements; ++elem) {
			DUMP("call build LS RHS. passing elem " << elem, DBG_DETAIL);
			auto &indA = elementToAtoms[elem];
			auto &indN = elementToNodes[elem];
			auto &mat = SubMatrix[elem];
			for (UInt i = 0; i < mat.rows() - 1; ++i) // loop on atoms
			for (UInt j = 0; j < mat.cols(); ++j) // loop on nodes
			for (UInt k = 0; k < stride; ++k) { // loop on dimensions
			rhs(indN[j], k) += mat(i, j) * v(indA[i], k);
			}
			}
			}
			/* --------------------------------------------------------------------------
			*/
			inline void ShapeMatrix::buildFormeFaibleCMatrix(Eigen::MatrixXf &C,
			std::vector<Real> &lambdas) {

			/* le produit de deux shape sur la position d'un atome
			n'est nul que si les deux shapes font partie
			d'un element commun : on a une seul boucle pricipale sur
			les elements */

			for (UInt k = 0; k < nb_elements; ++k) {
			auto &indN = elementToNodes[k];
			auto &indA = elementToAtoms[k];

			auto &mat = SubMatrix[k];

			/* ici on met la Real boucle sur les noeuds pour calculer
			la sous matrice de la grande */

			for (UInt j1 = 0; j1 < mat.cols(); ++j1)
			for (UInt j2 = 0; j2 < mat.cols(); ++j2) {
			/* ici on met la boucle de la somme sur les atomes */
			for (UInt i = 0; i < mat.rows() - 1; ++i) {

			C(indN[j1], indN[j2]) += mat(i, j1) * mat(i, j2) / lambdas[indA[i]];
			DUMP("C[" << indN[j1] << "," << indN[j2]
			<< "] = " << C(indN[j1], indN[j2]),
			DBG_ALL);
			}
			}
			}
			}
			/* --------------------------------------------------------------------------
			*/

			inline void ShapeMatrix::correctFormeFaible(
			[[gnu::unused]] Real &v, [[gnu::unused]] Array &rhs,
			[[gnu::unused]] std::vector<Real> &lamb, [[gnu::unused]] UInt at,
			[[gnu::unused]] UInt elem) {

			LM_TOIMPLEMENT;
			// UInt *indN = elementToNodes[elem];
			// UInt *indA = elementToAtoms[elem];
			// auto &mat = SubMatrix[elem];

			// /* trouve l'atome concerne */
			// UInt i, j;
			// for (i = 0; i < mat.rows() - 1; ++i)
			// if (indA[i] == at)
			// break;

			// for (j = 0; j < mat.cols(); ++j) {
			// v -= mat(i, j) * rhs(indN[j]) / lamb[at];
			// }
			}
			/* ------------------------------------------------------------------------ */
			inline void ShapeMatrix::swapAtomIndirections() {
			LM_TOIMPLEMENT;
			// UInt **temp = currentElementToAtom;
			// currentElementToAtom = newIndirectionAtom;
			// newIndirectionAtom = temp;
			elementToAtoms.swap();
			}
			/* ------------------------------------------------------------------------ */
			inline void ShapeMatrix::averageOnElements(ContainerArray<Real> &DataA,
			ContainerArray<Real> &DataM) {

			DataM.setZero();
			// loop over elements
			for (UInt elem = 0; elem < nb_elements; ++elem) {
			auto &indA = elementToAtoms[elem];
			auto &mat = SubMatrix[elem];
			DataM.row(elem).setZero();
			// loop over atoms
			for (UInt i = 0; i < mat.rows() - 1; ++i) {
			DataM.row(elem) += DataA.row(indA[i]);
			}
			DataM.row(elem) /= (mat.rows() - 1);
			}
			}
			/* --------------------------------------------------------------------------
			*/

			inline void ShapeMatrix::alterMatrixIndexesForPBC(
			std::vector<std::pair<UInt, UInt>> &pbc_pairs) {

			for (UInt k = 0; k < nb_elements; ++k) {
			auto &indN = elementToNodes[k];
			auto &mat = SubMatrix[k];
			for (UInt j = 0; j < mat.cols(); ++j) {
			for (UInt p = 0; p < pbc_pairs.size(); p++) {
			UInt i1 = pbc_pairs[p].first;
			UInt i2 = pbc_pairs[p].second;
			if (indN[j] == i1)
			indN[j] = i2;
			}
			}
			}
			}

			/* --------------------------------------------------------------------------
			*/
			__END_LIBMULTISCALE__

			#endif /* __LIBMULTISCALE_SHAPE_MATRIX_HH__ */

libmultiscale/src/coupling/shape_matrix.hh3a940948c399master

shape_matrix.hhView Options

libmultiscale/src/coupling/shape_matrix.hh
3a940948c399master

shape_matrix.hh
View Options