Homec4science
Diffusion LibMultiScale (features/coupler_component) Browse

libmultiscale/src/coupling/shape_matrix.hh
0b9417d37655features/coupler_compon…

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 <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;
};
/* ------------------------------------------------------------------------ */
template <UInt Dim> class ShapeMatrix {
public:
ShapeMatrix(UInt t);
~ShapeMatrix();
inline void setSub(UInt i, UInt m, UInt n);
inline void allocate();
inline void print();
inline void changeAtomIndirection(UInt el, UInt old_ind_atom,
UInt new_ind_atom);
inline void fill(UInt atom, UInt el, std::vector<UInt> &node_global,
std::vector<UInt> &nodes, std::vector<Real> &values);
inline void removeAtom(UInt at, UInt elem);
template <typename _Vec_>
inline Vector<Dim> interpolate(UInt at, _Vec_ &field, UInt elem);
inline Real interpolate(ContainerArray<Real> &atom_data,
ContainerArray<Real> &node_data);
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);
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 applyContinuumConstraint(_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);
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 auto indirectionAtomTable();
inline void extractShapeMatrix(math::Matrix *A);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
//! 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;
};
/* -------------------------------------------------------------------------- */
template <UInt Dim> ShapeMatrix<Dim>::ShapeMatrix(UInt nb_elements) {
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);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> ShapeMatrix<Dim>::~ShapeMatrix() {}
/* -------------------------------------------------------------------------- */
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);
matrix_sizes[i].first = m;
matrix_sizes[i].second = n;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void ShapeMatrix<Dim>::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;
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void ShapeMatrix<Dim>::print() {
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;
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void ShapeMatrix<Dim>::changeAtomIndirection(UInt el, UInt old_ind_atom,
UInt new_ind_atom) {
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
<< ")");
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void ShapeMatrix<Dim>::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;
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void ShapeMatrix<Dim>::removeAtom(UInt at, UInt elem) {
// 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;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void ShapeMatrix<Dim>::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);
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void ShapeMatrix<Dim>::buildThermalStuff(std::vector<Real> &nodal_out,
std::vector<Real> &atom_shape) {
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]];
}
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void ShapeMatrix<Dim>::buildContinuConstraint(
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>
template <typename _Vec_>
void ShapeMatrix<Dim>::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>
template <typename _Vec_>
void ShapeMatrix<Dim>::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>
template <typename _Vec_>
void ShapeMatrix<Dim>::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>
template <typename _Vec_>
Vector<Dim> ShapeMatrix<Dim>::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;
}
/* --------------------------------------------------------------------------
*/
template <UInt Dim>
Real ShapeMatrix<Dim>::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;
}
/* --------------------------------------------------------------------------
*/
template <UInt Dim>
void ShapeMatrix<Dim>::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] != -1)
atoms.push_back(indA[i]);
}
/* --------------------------------------------------------------------------
*/
template <UInt Dim>
void ShapeMatrix<Dim>::getRelatedNodesToElem(UInt elem,
std::vector<UInt> &nodes) {
auto &indGN = elementToNodesGlobalIndex[elem];
auto &mat = SubMatrix[elem];
nodes = indGN;
}
/* --------------------------------------------------------------------------
*/
template <UInt Dim> void ShapeMatrix<Dim>::extractShapeMatrix(math::Matrix *A) {
A->zero();
// 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);
}
}
}
}
/* --------------------------------------------------------------------------
*/
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 < 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");
A->printToFile(f);
}
/* --------------------------------------------------------------------------
*/
template <UInt Dim>
void ShapeMatrix<Dim>::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);
}
}
}
/* --------------------------------------------------------------------------
*/
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 < 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);
}
}
}
}
/* --------------------------------------------------------------------------
*/
template <UInt Dim>
void ShapeMatrix<Dim>::correctFormeFaible(Real &v, Array &rhs,
std::vector<Real> &lamb, UInt at,
UInt elem) {
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];
}
}
/* ------------------------------------------------------------------------ */
template <UInt Dim> void ShapeMatrix<Dim>::swapAtomIndirections() {
// UInt **temp = currentElementToAtom;
// currentElementToAtom = newIndirectionAtom;
// newIndirectionAtom = temp;
elementToAtoms.swap();
}
/* ------------------------------------------------------------------------ */
template <UInt Dim>
void ShapeMatrix<Dim>::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);
}
}
/* --------------------------------------------------------------------------
*/
template <UInt Dim>
void ShapeMatrix<Dim>::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__ */
c4science · Help