Homec4science
Diffusion LibMultiScale (features/pybind) Browse

libmultiscale/src/coupling/shape_matrix.hh
b05d5f2aba78features/pybind

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 "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(std::vector<Real> & DataA,std::vector<Real> & DataM,
UInt stride);
inline void alterMatrixIndexesForPBC(std::vector<std::pair<UInt, UInt> > & pbc_pairs);
/* ------------------------------------------------------------------------ */
/* Methods (specific to a coupling type) */
/* ------------------------------------------------------------------------ */
inline void buildNodeShape(std::vector<Real> & node_shape);
inline void buildThermalStuff(std::vector<Real> & nodal_out,
std::vector<Real> & atom_shape);
inline void buildContinuConstraint(std::vector<Real> & A,std::vector<Real> & node_shape,
std::vector<Real> & node_weight);
template <typename _Vec_>
inline void buildRHS(std::vector<Real> & rhs,_Vec_ & v);
template <typename _Vec_>
inline void applyChanging(_Vec_ & v,std::vector<Real> & multL,std::vector<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(std::vector<Real> & rhs,std::vector<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(std::vector<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(std::vector<Real> & A,
std::vector<Real> & node_shape,
std::vector<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(std::vector<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,std::vector<Real> & multL,
std::vector<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(std::vector<Real> & rhs,
std::vector<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(std::vector<Real> & DataA,
std::vector<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