Homec4science
Diffusion LibMultiScale (features/reshape_compute_python) Browse

libmultiscale/src/coupling/shape_matrix.hh
27ca80ef7f29features/reshape_comput…

Last Change
Raw File

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__ */
c4science · Help