shape_matrix.hh
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Wed, Nov 13, 01:29

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 "full_matrix.hh"
	/* -------------------------------------------------------------------------- */

	__BEGIN_LIBMULTISCALE__

	template <UInt Dim> class ShapeMatrix {
	public:
	/* ------------------------------------------------------------------------ */
	/* Constructors/Destructors */
	/* ------------------------------------------------------------------------ */

	ShapeMatrix(UInt t);
	~ShapeMatrix();

	/* ------------------------------------------------------------------------ */
	/* Methods (management) */
	/* ------------------------------------------------------------------------ */

	inline void setSub(UInt i, UInt m, UInt n);
	inline void allocate();
	inline void print();
	inline void changeAtomIndirection(UInt el, UInt old_ind_atome,
	UInt new_ind_atome);
	inline void fill(UInt atome, UInt node_global, UInt node, UInt el,
	Real value);
	inline void removeAtom(UInt at, UInt elem);
	template <typename _Vec_>
	inline Real interpolate(UInt at, _Vec_ &field, UInt elem, UInt deg);
	inline Real interpolate(std::vector<Real> &atom_data,
	std::vector<Real> &node_data, UInt stride);
	inline void getRelatedAtomsToElem(UInt elem, std::vector<UInt> &atoms);
	inline void getRelatedNodesToElem(UInt elem, std::vector<UInt> &nodes);
	inline void swapAtomIndirections();
	inline void averageOnElements(ContainerArray<Real> &DataA,
	ContainerArray<Real> &DataM, UInt stride);

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

	/* ------------------------------------------------------------------------ */
	/* 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 buildContinuConstraint(ContainerArray<Real> &A,
	ContainerArray<Real> &node_shape,
	ContainerArray<Real> &node_weight);
	template <typename _Vec_>
	inline void buildRHS(ContainerArray<Real> &rhs, _Vec_ &v);
	template <typename _Vec_>
	inline void applyChanging(_Vec_ &v, ContainerArray<Real> &multL,
	ContainerArray<Real> &lambdas);
	template <typename _Vec_>
	inline void applyChangingThermal(_Vec_ &vDisp, _Vec_ &vVel, _Vec_ &vAcc,
	_Vec_ &masses, std::vector<Real> &multL,
	Real deltat);
	inline void buildLeastSquareMatrix(math::Matrix *A);
	inline void buildLeastSquareRHS(ContainerArray<Real> &rhs,
	ContainerArray<Real> &v, UInt stride);
	inline void buildFormeFaibleCMatrix(math::Matrix &C,
	std::vector<Real> &lambdas);
	inline void correctFormeFaible(Real &v, Array &rhs, std::vector<Real> &lamb,
	UInt at, UInt elem);

	/* ------------------------------------------------------------------------ */
	/* Accessors */
	/* ------------------------------------------------------------------------ */

	inline UInt **indirectionAtomTable();
	inline void extractShapeMatrix(math::Matrix *A);

	/* ------------------------------------------------------------------------ */
	/* Class Members */
	/* ------------------------------------------------------------------------ */

	private:
	//! number of stored elements
	UInt taille;
	//! size of submatrices (m1,n1,m2,n2,...)
	UInt *tailles;
	//! array of submatrices
	math::Matrix **SubMatrix;
	//! index map to atoms (current)
	UInt **currentIndirectionAtome;
	//! index map to atoms (second : swapped system)
	UInt **newIndirectionAtome;
	//! index map to atoms
	UInt **indirectionAtome[2];
	//! index map to nodes
	UInt **indirectionNode;
	//! index map to node in the global index (big vector)
	UInt **indirectionGlobalNode;
	//! indicator of filling for \phi_I(X_i) submatrices for atoms
	UInt *remplissageA;
	//! indicator of filling for \phi_I(X_i) submatrices for nodes
	UInt *remplissageN;
	//! contiguous array of data where to store submatrices
	Real *tab_contigu;
	//! index map to atomic arrays (contiguous)
	UInt *indA_contigu[2];
	//! index map to nodal arrays (contiguous)
	UInt *indN_contigu;
	//! global index map to nodal arrays (contiguous)
	UInt *indGN_contigu;
	};

	/* -------------------------------------------------------------------------- */
	template <UInt Dim> ShapeMatrix<Dim>::ShapeMatrix(UInt t) {
	DUMP("creating shapematrix for " << t << " elements", DBG_WARNING);
	SubMatrix = new math::Matrix *[t];
	tailles = new UInt[2 * t];
	indirectionAtome[0] = new UInt *[t];
	indirectionAtome[1] = new UInt *[t];
	memset(indirectionAtome[0], 0, sizeof(int ) t);
	memset(indirectionAtome[1], 0, sizeof(int ) t);
	currentIndirectionAtome = indirectionAtome[0];
	newIndirectionAtome = indirectionAtome[1];
	indirectionNode = new UInt *[t];
	indirectionGlobalNode = new UInt *[t];
	remplissageA = new UInt[t];
	remplissageN = new UInt[t];
	memset(remplissageA, 0, sizeof(UInt) * t);
	memset(remplissageN, 0, sizeof(UInt) * t);

	tab_contigu = NULL;
	taille = t;
	};
	/* -------------------------------------------------------------------------- */
	template <UInt Dim> ShapeMatrix<Dim>::~ShapeMatrix() {
	for (UInt i = 0; i < taille; ++i)
	delete SubMatrix[i];

	if (tab_contigu)
	delete[] tab_contigu;
	if (indA_contigu[0])
	delete[] indA_contigu[0];
	if (indA_contigu[1])
	delete[] indA_contigu[1];
	if (indN_contigu)
	delete[] indN_contigu;
	if (indGN_contigu)
	delete[] indGN_contigu;

	delete[] remplissageA;
	delete[] remplissageN;
	delete[] indirectionNode;
	delete[] indirectionGlobalNode;
	delete[] indirectionAtome[0];
	delete[] indirectionAtome[1];
	delete[] SubMatrix;
	delete[] tailles;
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim> void ShapeMatrix<Dim>::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);
	tailles[i * 2] = m;
	tailles[i * 2 + 1] = n;
	};
	/* -------------------------------------------------------------------------- */
	template <UInt Dim> void ShapeMatrix<Dim>::allocate() {
	// premier passage pour calculer tout les
	// coeficients de la smatrice globale

	UInt m, n;

	UInt total_size = 0;
	UInt total_atoms = 0;
	UInt total_nodes = 0;

	for (UInt i = 0; i < taille; ++i) {
	m = tailles[i * 2];
	n = tailles[i * 2 + 1];

	total_size += (m + 1) * n;
	total_atoms += m;
	total_nodes += n;
	}

	// allocation
	tab_contigu = new Real[total_size];
	UInt index = 0;

	indA_contigu[0] = new UInt[total_atoms];
	indA_contigu[1] = new UInt[total_atoms];
	for (UInt i = 0; i < total_atoms; ++i) {
	indA_contigu[0][i] = UINT_MAX;
	indA_contigu[1][i] = UINT_MAX;
	}
	UInt index_atoms = 0;

	indN_contigu = new UInt[total_nodes];
	indGN_contigu = new UInt[total_nodes];
	UInt index_nodes = 0;

	DUMP("ShapeMatrix allocate "
	<< sizeof(Real) * total_size +
	2 * (total_atoms + total_nodes) * sizeof(UInt)
	<< " bytes",
	DBG_INFO_STARTUP);

	// deuxieme passage pour allouer les structures de control des matrices
	for (UInt i = 0; i < taille; ++i) {
	m = tailles[i * 2];
	n = tailles[i * 2 + 1];

	SubMatrix[i] = new math::Matrix(m + 1, n);

	indirectionAtome[0][i] = indA_contigu[0] + index_atoms;
	indirectionAtome[1][i] = indA_contigu[1] + index_atoms;
	indirectionNode[i] = indN_contigu + index_nodes;
	indirectionGlobalNode[i] = indGN_contigu + index_nodes;

	// incrementation des indexs
	index += (m + 1) * n;
	index_atoms += m;
	index_nodes += n;
	}
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim> void ShapeMatrix<Dim>::print() {
	for (UInt k = 0; k < taille; ++k) {
	if (!remplissageA[k])
	return;
	math::Matrix &mat = *SubMatrix[k];

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

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

	std::cout << std::endl;

	SubMatrix[k]->print();
	}
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim>
	void ShapeMatrix<Dim>::changeAtomIndirection(UInt el, UInt old_ind_atome,
	UInt new_ind_atome) {
	math::Matrix &mat = *SubMatrix[el];
	UInt *indA = currentIndirectionAtome[el];
	UInt *indAN = newIndirectionAtome[el];

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

	UInt i = 0;
	for (i = 0; i < mat.m() - 1; ++i) {
	if (indA[i] == old_ind_atome) {
	indAN[i] = new_ind_atome;
	break;
	}
	}
	LM_ASSERT(i != mat.m() - 1,
	"atom not found this should not append (index was " << old_ind_atome
	<< ")");
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim>
	void ShapeMatrix<Dim>::fill(UInt atome, UInt node_global, UInt node, UInt el,
	Real value) {

	LM_ASSERT(!std::isnan(value) && !std::isinf(value),
	"the passed value for atom " << atome << " node_global "
	<< node_global << " node " << node
	<< "el " << el);

	UInt *indA = currentIndirectionAtome[el];
	indA[remplissageA[el]] = atome;

	UInt *indN = indirectionNode[el];
	indN[remplissageN[el]] = node;

	DUMP("setting shapes for atome " << atome << " and element " << el << " node "
	<< node << " node_global " << node_global
	<< " value " << value,
	DBG_ALL);
	UInt *indGN = indirectionGlobalNode[el];
	indGN[remplissageN[el]] = node_global;
	// je prend la matrice dense qui correspond a mon element
	math::Matrix &m = *SubMatrix[el];
	m(remplissageA[el], remplissageN[el]) = value;
	// fait les sommes partielles sur les noeuds
	m(m.m() - 1, remplissageN[el]) += value;

	++remplissageN[el];
	if (remplissageN[el] == m.n()) {
	remplissageN[el] = 0;
	++remplissageA[el];
	}
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim> void ShapeMatrix<Dim>::removeAtom(UInt at, UInt elem) {
	// I don't really remove the atom.
	// I only substract the contribution to the cumulative part at bottom of
	// submatrix.
	// since assoc will be modified it is sufficient (I think !)
	math::Matrix &mat = *SubMatrix[elem];
	UInt *indA = currentIndirectionAtome[elem];
	UInt i = 0;
	for (i = 0; i < mat.m() - 1; ++i)
	if (indA[i] == at)
	break;

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

	for (UInt j = 0; j < mat.n(); ++j) {
	mat(mat.m() - 1, j) -= mat(i, j);
	}
	indA[i] = UINT_MAX;
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim>
	void ShapeMatrix<Dim>::buildNodeShape(ContainerArray<Real> &node_shape) {
	for (UInt k = 0; k < taille; ++k) {
	UInt *indN = indirectionNode[k];
	math::Matrix &mat = *SubMatrix[k];
	for (UInt j = 0; j < mat.n(); ++j) {
	DUMP("Summing node_shape[" << indN[j] << " (el= " << k << ", j = " << j
	<< ")] with " << mat(mat.m() - 1, j),
	DBG_ALL);
	node_shape[indN[j]] += mat(mat.m() - 1, j);
	}
	}
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim>
	void ShapeMatrix<Dim>::buildThermalStuff(std::vector<Real> &nodal_out,
	std::vector<Real> &atom_shape) {
	for (UInt elem = 0; elem < taille; ++elem) {
	UInt *indN = indirectionNode[elem];
	UInt *indA = currentIndirectionAtome[elem];
	math::Matrix &mat = *SubMatrix[elem];
	for (UInt i = 0; i < mat.m() - 1; ++i) {
	for (UInt j = 0; j < mat.n(); ++j) {
	nodal_out[indN[j]] += mat(i, j) * atom_shape[indA[i]];
	}
	}
	}
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim>
	void ShapeMatrix<Dim>::buildContinuConstraint(
	ContainerArray<Real> &A, ContainerArray<Real> &node_shape,
	ContainerArray<Real> &node_weight) {
	for (UInt k = 0; k < taille; ++k) {
	UInt *indA = currentIndirectionAtome[k];
	UInt *indN = indirectionNode[k];
	math::Matrix &mat = *SubMatrix[k];
	for (UInt i = 0; i < mat.m() - 1; ++i)
	if (indA[i] != UINT_MAX)
	for (UInt j = 0; j < mat.n(); ++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>
	template <typename _Vec_>
	void ShapeMatrix<Dim>::buildRHS(ContainerArray<Real> &rhs, _Vec_ &v) {
	// je parcours tout les blocs
	for (UInt k = 0; k < taille; ++k) {
	// je chope la sous matrice
	math::Matrix &mat = *SubMatrix[k];
	// je chope les indirection sur les noeuds globaux comme
	// j'ai un acces direct au grand vecteur
	UInt *indGN = indirectionGlobalNode[k];
	UInt *indA = currentIndirectionAtome[k];
	// je parcours tout les atomes du bloc
	for (UInt i = 0; i < mat.m() - 1; ++i)
	if (indA[i] != UINT_MAX)
	// je parcours tout les noeuds
	for (UInt j = 0; j < mat.n(); ++j) {
	DUMP("Before rhs[" << indA[i] << "]=" << rhs[Dim * indA[i] + 1],
	DBG_ALL);
	rhs[Dim * indA[i]] += v[indGN[j] * Dim] * mat(i, j);
	if (Dim > 1)
	rhs[Dim * indA[i] + 1] += v[indGN[j] * Dim + 1] * mat(i, j);
	if (Dim == 3)
	rhs[Dim * indA[i] + 2] += v[indGN[j] * Dim + 2] * mat(i, j);
	DUMP("build rhs to " << Dim * indA[i] << " " << Dim * indA[i] + 1
	<< " " << Dim * indA[i] + 3,
	DBG_ALL);
	DUMP("Building " << rhs[Dim * indA[i] + 1] << "=rhs[" << indA[i]
	<< "] with " << mat(i, j) << "*"
	<< v[indGN[j] * Dim + 1] << " = "
	<< v[indGN[j] * Dim + 1] * mat(i, j)
	<< " (index = " << indGN[j] << ")",
	DBG_ALL);
	}
	}
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim>
	template <typename _Vec_>
	void ShapeMatrix<Dim>::applyChanging(_Vec_ &v, ContainerArray<Real> &multL,
	ContainerArray<Real> &lambdas) {
	// je parcours tout les blocs
	for (UInt k = 0; k < taille; ++k) {
	// je chope la sous matrice
	math::Matrix &mat = *SubMatrix[k];
	// je chope les indirection sur les noeuds globaux comme
	// j'ai un acces direct au grand vecteur
	UInt *indGN = indirectionGlobalNode[k];
	UInt *indN = indirectionNode[k];
	UInt *indA = currentIndirectionAtome[k];
	// je parcours tout les atomes du bloc
	for (UInt i = 0; i < mat.m() - 1; ++i)
	if (indA[i] != UINT_MAX)
	// je parcours tout les noeuds
	for (UInt j = 0; j < mat.n(); ++j) {
	DUMP("adding v[" << indGN[j] << "] with " << -1.0 * mat(i, j) << " * "
	<< multL[indA[i] * Dim + 1] << " {rhs[" << indA[i]
	<< "]}/ " << lambdas[indN[j]] << " = "
	<< -1.0 * mat(i, j) * multL[indA[i] * Dim + 1] /
	lambdas[indN[j]],
	DBG_ALL);

	v.add(indGN[j] * Dim,
	-1.0 * mat(i, j) * multL[indA[i] * Dim] / lambdas[indN[j]]);
	if (Dim > 1)
	v.add(indGN[j] * Dim + 1, -1.0 * mat(i, j) *
	multL[indA[i] * Dim + 1] /
	lambdas[indN[j]]);
	if (Dim == 3)
	v.add(indGN[j] * Dim + 2, -1.0 * mat(i, j) *
	multL[indA[i] * Dim + 2] /
	lambdas[indN[j]]);
	}
	}
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim>
	template <typename _Vec_>
	void ShapeMatrix<Dim>::applyChangingThermal(_Vec_ &vDisp, _Vec_ &vVel,
	_Vec_ &vAcc, _Vec_ &masses,
	std::vector<Real> &multL,
	Real deltat) {
	// je parcours tout les blocs
	for (UInt k = 0; k < taille; ++k) {
	// je chope la sous matrice
	math::Matrix &mat = *SubMatrix[k];
	// je chope les indirection sur les noeuds globaux comme
	// j'ai un acces direct au grand vecteur
	UInt *indGN = indirectionGlobalNode[k];
	// UInt * indN = indirectionNode[k];
	UInt *indA = currentIndirectionAtome[k];
	// je parcours tout les atomes du bloc
	for (UInt i = 0; i < mat.m() - 1; ++i)
	if (indA[i] != UINT_MAX)
	// je parcours tout les noeuds
	for (UInt j = 0; j < mat.n(); ++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>
	template <typename _Vec_>
	Real ShapeMatrix<Dim>::interpolate(UInt at, _Vec_ &field, UInt elem, UInt deg) {
	Real res = 0;
	LM_ASSERT(elem < taille, "overflow detected");
	UInt *indA = currentIndirectionAtome[elem];
	UInt *indGN = indirectionGlobalNode[elem];
	math::Matrix &mat = *SubMatrix[elem];

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

	for (UInt elem = 0; elem < taille; ++elem) {
	UInt *indN = indirectionNode[elem];
	UInt *indA = currentIndirectionAtome[elem];
	math::Matrix &mat = *SubMatrix[elem];
	for (UInt i = 0; i < mat.m() - 1; ++i) {
	for (UInt j = 0; j < mat.n(); ++j) {
	for (UInt k = 0; k < stride; ++k) {
	atom_data[stride * indA[i] + k] +=
	mat(i, j) * node_data[indN[j] * stride + k];
	}
	}
	}
	}
	return 0;
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim>
	void ShapeMatrix<Dim>::getRelatedAtomsToElem(UInt elem,
	std::vector<UInt> &atoms) {

	UInt *indA = currentIndirectionAtome[elem];
	math::Matrix &mat = *SubMatrix[elem];

	atoms.clear();

	// parcours les atomes
	UInt i;
	for (i = 0; i < mat.m() - 1; ++i)
	if (indA[i] != -1)
	atoms.push_back(indA[i]);
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim>
	void ShapeMatrix<Dim>::getRelatedNodesToElem(UInt elem,
	std::vector<UInt> &nodes) {

	UInt *indGN = indirectionGlobalNode[elem];
	math::Matrix &mat = *SubMatrix[elem];

	nodes.clear();

	// parcours les noeuds
	UInt j;
	for (j = 0; j < mat.n(); ++j)
	nodes.push_back(indGN[j]);
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim> void ShapeMatrix<Dim>::extractShapeMatrix(math::Matrix *A) {
	A->zero();
	// I loop over elements
	for (UInt k = 0; k < taille; ++k) {
	UInt *indA = currentIndirectionAtome[k];
	UInt *indN = indirectionNode[k];
	math::Matrix &mat = *SubMatrix[k];
	// loop over nodes
	for (UInt j = 0; j < mat.n(); ++j) {
	// loop over atoms
	for (UInt i = 0; i < mat.m() - 1; ++i) {
	(*A)(indA[i], indN[j]) += mat(i, j);
	}
	}
	}
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim>
	void ShapeMatrix<Dim>::buildLeastSquareMatrix(math::Matrix *A) {
	A->zero();

	/* 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 < taille; ++k) {
	UInt *indN = indirectionNode[k];
	math::Matrix &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.n(); ++j1)
	for (UInt j2 = 0; j2 < mat.n(); ++j2) {
	/* ici on met la boucle de la somme sur les atomes */
	for (UInt i = 0; i < mat.m() - 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");
	A->printToFile(f);
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim>
	void ShapeMatrix<Dim>::buildLeastSquareRHS(ContainerArray<Real> &rhs,
	ContainerArray<Real> &v,
	UInt stride) {

	DUMP("stride is " << stride << " taille " << taille, DBG_DETAIL);
	// loop on elements
	for (UInt elem = 0; elem < taille; ++elem) {
	DUMP("call build LS RHS. passing elem " << elem, DBG_DETAIL);
	UInt *indA = currentIndirectionAtome[elem];
	UInt *indN = indirectionNode[elem];
	math::Matrix &mat = *SubMatrix[elem];
	for (UInt i = 0; i < mat.m() - 1; ++i) // loop on atoms
	for (UInt j = 0; j < mat.n(); ++j) // loop on nodes
	for (UInt k = 0; k < stride; ++k) { // loop on dimensions
	rhs[indN[j] * stride + k] += mat(i, j) * v[indA[i] * stride + k];
	}
	}
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim>
	void ShapeMatrix<Dim>::buildFormeFaibleCMatrix(math::Matrix &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 < taille; ++k) {
	UInt *indN = indirectionNode[k];
	UInt *indA = currentIndirectionAtome[k];

	math::Matrix &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.n(); ++j1)
	for (UInt j2 = 0; j2 < mat.n(); ++j2) {
	/* ici on met la boucle de la somme sur les atomes */
	for (UInt i = 0; i < mat.m() - 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);
	}
	}
	}
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim>
	void ShapeMatrix<Dim>::correctFormeFaible(Real &v, Array &rhs,
	std::vector<Real> &lamb, UInt at,
	UInt elem) {

	UInt *indN = indirectionNode[elem];
	UInt *indA = currentIndirectionAtome[elem];
	math::Matrix &mat = *SubMatrix[elem];

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

	for (j = 0; j < mat.n(); ++j) {
	v -= mat(i, j) * rhs(indN[j]) / lamb[at];
	}
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim> UInt **ShapeMatrix<Dim>::indirectionAtomTable() {
	return indA_contigu;
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim> void ShapeMatrix<Dim>::swapAtomIndirections() {
	UInt **temp = currentIndirectionAtome;
	currentIndirectionAtome = newIndirectionAtome;
	newIndirectionAtome = temp;
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim>
	void ShapeMatrix<Dim>::averageOnElements(ContainerArray<Real> &DataA,
	ContainerArray<Real> &DataM,
	UInt stride) {
	// I loop over elements
	for (UInt elem = 0; elem < taille; ++elem) {
	UInt *indA = currentIndirectionAtome[elem];
	// UInt * indN = indirectionNode[elem];
	math::Matrix &mat = *SubMatrix[elem];
	DataM[elem * stride] = 0;
	DataM[elem * stride + 1] = 0;
	DataM[elem * stride + 2] = 0;
	// loop over atoms
	for (UInt i = 0; i < mat.m() - 1; ++i) {
	for (UInt k = 0; k < stride; ++k) {
	DataM[elem * stride + k] += DataA[indA[i] * stride + k];
	DUMP("adding " << DataA[indA[i] * stride + k] << " to elem " << elem
	<< "(" << k << ")",
	DBG_ALL);
	}
	}
	for (UInt k = 0; k < stride; ++k)
	DataM[elem * stride + k] /= (mat.m() - 1);
	}
	}
	/* -------------------------------------------------------------------------- */

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

	for (UInt k = 0; k < taille; ++k) {
	UInt *indN = indirectionNode[k];
	math::Matrix &mat = *SubMatrix[k];
	for (UInt j = 0; j < mat.n(); ++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__ */

shape_matrix.hhNo OneTemporaryActions

File Metadata

shape_matrix.hhView Options

Event Timeline

shape_matrix.hh
No OneTemporary
Actions

shape_matrix.hh
View Options