Page MenuHomec4science
Files F67870971

shape_matrix.hh
No OneTemporary
Actions

File Metadata

Created
Mon, Jun 24, 19:53

shape_matrix.hh
View Options

/**
* @file shape_matrix.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Fri Jul 11 15:47:44 2014
*
* @brief Data structure for the storage of shape function values on
* atomic/point positions
*
* @section LICENSE
*
* Copyright INRIA and CEA
*
* The LibMultiScale is a C++ parallel framework for the multiscale
* coupling methods dedicated to material simulations. This framework
* provides an API which makes it possible to program coupled simulations
* and integration of already existing codes.
*
* This Project was initiated in a collaboration between INRIA Futurs Bordeaux
* within ScAlApplix team and CEA/DPTA Ile de France.
* The project is now continued at the Ecole Polytechnique Fédérale de Lausanne
* within the LSMS/ENAC laboratory.
*
* This software is governed by the CeCILL-C license under French law and
* abiding by the rules of distribution of free software. You can use,
* modify and/ or redistribute the software under the terms of the CeCILL-C
* license as circulated by CEA, CNRS and INRIA at the following URL
* "http://www.cecill.info".
*
* As a counterpart to the access to the source code and rights to copy,
* modify and redistribute granted by the license, users are provided only
* with a limited warranty and the software's author, the holder of the
* economic rights, and the successive licensors have only limited
* liability.
*
* In this respect, the user's attention is drawn to the risks associated
* with loading, using, modifying and/or developing or reproducing the
* software by the user in light of its specific status of free software,
* that may mean that it is complicated to manipulate, and that also
* therefore means that it is reserved for developers and experienced
* professionals having in-depth computer knowledge. Users are therefore
* encouraged to load and test the software's suitability as regards their
* requirements in conditions enabling the security of their systems and/or
* data to be ensured and, more generally, to use and operate it in the
* same conditions as regards security.
*
* The fact that you are presently reading this means that you have had
* knowledge of the CeCILL-C license and that you accept its terms.
*
*/
#ifndef __LIBMULTISCALE_SHAPE_MATRIX_HH__
#define __LIBMULTISCALE_SHAPE_MATRIX_HH__
/* -------------------------------------------------------------------------- */
#include "container_array.hh"
/* ------------------------------------------------------------------------ */
#include <Eigen/Sparse>
/* ------------------------------------------------------------------------ */
__BEGIN_LIBMULTISCALE__
/* ------------------------------------------------------------------------ */
template <typename T> struct SwapableObject {
SwapableObject() {
current = &obj1;
_next = &obj2;
}
void swap() {
auto *temp = current;
current = _next;
_next = temp;
}
T &operator*() { return *current; };
T *operator->() { return current; };
template <typename T1> decltype(auto) operator[](T1 &&idx) {
return (*current)[idx];
}
T &next() { return *_next; }
T obj1;
T obj2;
T *current;
T *_next;
};
/* ------------------------------------------------------------------------ */
class ShapeMatrix
: public Eigen::SparseMatrix<libmultiscale::Real, Eigen::RowMajor> {
public:
using Eigen::SparseMatrix<libmultiscale::Real, Eigen::RowMajor>::SparseMatrix;
~ShapeMatrix();
inline Real interpolate(ContainerArray<Real> &atom_data,
ContainerArray<Real> &node_data);
inline void removeAtom(UInt at, UInt elem);
inline void changeAtomIndirection(UInt el, UInt old_ind_atom,
UInt new_ind_atom);
inline void swapAtomIndirections();
//! return true if the matrix is ready to be used
bool isFilled();
private:
// inline void setSub(UInt i, UInt m, UInt n);
// inline void allocate();
// inline void print();
// inline void fill(UInt atom, UInt el, std::vector<UInt> &node_global,
// std::vector<UInt> &nodes, std::vector<Real> &values);
// template <UInt Dim, typename _Vec_>
// inline Vector<Dim> interpolate(UInt at, _Vec_ &field, UInt elem);
inline void getRelatedAtomsToElem(UInt elem, std::vector<UInt> &atoms);
inline void getRelatedNodesToElem(UInt elem, std::vector<UInt> &nodes);
inline void averageOnElements(ContainerArray<Real> &DataA,
ContainerArray<Real> &DataM);
inline void
alterMatrixIndexesForPBC(std::vector<std::pair<UInt, UInt>> &pbc_pairs);
public:
/* ------------------------------------------------------------------------ */
/* Methods (specific to a coupling type) */
/* ------------------------------------------------------------------------ */
inline void buildNodeShape(ContainerArray<Real> &node_shape);
inline void buildThermalStuff(std::vector<Real> &nodal_out,
std::vector<Real> &atom_shape);
inline void buildContinuumConstraint(ContainerArray<Real> &A,
ContainerArray<Real> &node_shape,
ContainerArray<Real> &node_weight);
template <UInt Dim, typename Vec>
void buildRHS(ContainerArray<Real> &rhs, Vec &v);
template <UInt Dim, typename _Vec_>
inline void applyContinuumConstraint(_Vec_ &v, ContainerArray<Real> &multL,
ContainerArray<Real> &lambdas);
template <UInt Dim, typename _Vec_>
inline void applyChangingThermal(_Vec_ &vDisp, _Vec_ &vVel, _Vec_ &vAcc,
_Vec_ &masses, std::vector<Real> &multL,
Real deltat);
inline void buildLeastSquareMatrix(Eigen::MatrixXf &A);
inline void buildLeastSquareRHS(ContainerArray<Real> &rhs,
ContainerArray<Real> &v);
inline void buildFormeFaibleCMatrix(Eigen::MatrixXf &C,
std::vector<Real> &lambdas);
inline void correctFormeFaible(Real &v, Array &rhs, std::vector<Real> &lamb,
UInt at, UInt elem);
inline auto indirectionAtomTable();
inline void extractShapeMatrix(Eigen::MatrixXf &A);
private:
//! filled indicator
bool is_filled = false;
//! number of nodes in the mesh
UInt nb_nodes;
//! number of points
UInt nb_points;
//! number of stored elements
UInt nb_elements;
//! size of submatrices (m1,n1,m2,n2,...)
std::vector<std::pair<UInt, UInt>> matrix_sizes;
//! array of submatrices
std::vector<Eigen::MatrixXd> SubMatrix;
//! map element index to contained atom indexes
SwapableObject<std::vector<std::vector<UInt>>> elementToAtoms;
//! map element index to contained node indexes
std::vector<std::vector<UInt>> elementToNodes;
//! index map to node in the global index (big vector)
std::vector<std::vector<UInt>> elementToNodesGlobalIndex;
};
/* -------------------------------------------------------------------------- */
// inline ShapeMatrix::ShapeMatrix(UInt nb_elements, UInt nb_nodes,
// UInt nb_points) {
// this->nb_elements = nb_elements;
// this->nb_points = nb_points;
// this->nb_nodes = nb_nodes;
// DUMP("creating shapematrix for " << nb_elements << " elements",
// DBG_WARNING); SubMatrix.resize(nb_elements);
// matrix_sizes.resize(nb_elements);
// elementToAtoms.obj1.resize(nb_elements);
// elementToAtoms.obj2.resize(nb_elements);
// elementToNodes.resize(nb_elements);
// elementToNodesGlobalIndex.resize(nb_elements);
// }
/* -------------------------------------------------------------------------- */
// inline ShapeMatrix::ShapeMatrix(UInt nb_nodes, UInt nb_points) {
// this->nb_points = nb_points;
// this->nb_nodes = nb_nodes;
// DUMP("creating shapematrix for " << nb_elements << " elements",
// DBG_WARNING); SubMatrix.resize(nb_elements);
// matrix_sizes.resize(nb_elements);
// elementToAtoms.obj1.resize(nb_elements);
// elementToAtoms.obj2.resize(nb_elements);
// elementToNodes.resize(nb_elements);
// elementToNodesGlobalIndex.resize(nb_elements);
// }
// /* --------------------------------------------------------------------------
// */
// inline ShapeMatrix::ShapeMatrix(UInt nb_elements) {
// LM_TOIMPLEMENT;
// this->nb_elements = nb_elements;
// DUMP("creating shapematrix for " << nb_elements << " elements",
// DBG_WARNING); SubMatrix.resize(nb_elements);
// matrix_sizes.resize(nb_elements);
// elementToAtoms.obj1.resize(nb_elements);
// elementToAtoms.obj2.resize(nb_elements);
// elementToNodes.resize(nb_elements);
// elementToNodesGlobalIndex.resize(nb_elements);
// }
/* -------------------------------------------------------------------------- */
inline ShapeMatrix::~ShapeMatrix() {}
/* -------------------------------------------------------------------------- */
inline bool ShapeMatrix::isFilled() { return this->is_filled; }
/* -------------------------------------------------------------------------- */
// inline void ShapeMatrix::setSub(UInt i, UInt m, UInt n) {
// LM_ASSERT(m > 0 && n > 0,
// "can not set a smatrix bloc size of zero : m=" << m << " n=" <<
// n);
// matrix_sizes[i].first = m;
// matrix_sizes[i].second = n;
// }
/* -------------------------------------------------------------------------- */
// inline void ShapeMatrix::allocate() {
// // first pass to compute the global ShapeMatrix coefficients
// UInt m, n;
// UInt total_atoms = 0;
// UInt total_nodes = 0;
// for (UInt i = 0; i < nb_elements; ++i) {
// total_atoms += matrix_sizes[i].first;
// total_nodes += matrix_sizes[i].second;
// }
// UInt index_atoms = 0;
// UInt index_nodes = 0;
// // second pass to allocate the data structures of matrices
// for (UInt i = 0; i < nb_elements; ++i) {
// m = matrix_sizes[i].first;
// n = matrix_sizes[i].second;
// SubMatrix[i].resize(m + 1, n);
// SubMatrix[i].setZero();
// // increment of indexes
// index_atoms += m;
// index_nodes += n;
// }
// }
/* -------------------------------------------------------------------------- */
// inline void ShapeMatrix::print() {
// std::cout << *this << std::endl;
// for (UInt k = 0; k < nb_elements; ++k) {
// auto &mat = SubMatrix[k];
// if (mat.rows() == 1)
// return;
// DUMP("affiche les infos de " << k, DBG_DETAIL);
// for (UInt i = 0; i < mat.rows() - 1; ++i)
// std::cout << elementToAtoms[k][i] << "\t";
// std::cout << std::endl;
// std::cout << SubMatrix[k] << std::endl;
// }
//}
/* -------------------------------------------------------------------------- */
inline void ShapeMatrix::changeAtomIndirection(UInt el, UInt old_ind_atom,
UInt new_ind_atom) {
auto &mat = SubMatrix[el];
auto &indA = elementToAtoms[el];
auto &indAN = elementToAtoms.next()[el];
indAN.resize(indA.size());
// should be done in swapatomindirection TODO !!
for (UInt i = 0; i < mat.rows() - 1; ++i)
if (indA[i] == UINT_MAX)
indAN[i] = UINT_MAX;
UInt i = 0;
for (i = 0; i < mat.rows() - 1; ++i) {
if (indA[i] == old_ind_atom) {
indAN[i] = new_ind_atom;
break;
}
}
LM_ASSERT(i != mat.rows() - 1,
"atom not found this should not append (index was " << old_ind_atom
<< ")");
}
/* --------------------------------------------------------------------------
*/
// inline void ShapeMatrix::fill(UInt atom, UInt el,
// std::vector<UInt> &nodes_global,
// std::vector<UInt> &nodes,
// std::vector<Real> &values) {
// UInt atom_index = elementToAtoms[el].size();
// elementToAtoms[el].push_back(atom);
// if (elementToNodes[el].size() == 0) {
// elementToNodes[el] = nodes;
// elementToNodesGlobalIndex[el] = nodes_global;
// }
// for (UInt nd = 0; nd < nodes.size(); ++nd) {
// Real &value = values[nd];
// UInt &node_global [[gnu::unused]] = nodes_global[nd];
// UInt &node [[gnu::unused]] = nodes[nd];
// LM_ASSERT(!std::isnan(value) && !std::isinf(value),
// "the passed value for atom " << atom << " node_global "
// << node_global << " node " << node
// << "el " << el);
// DUMP("setting shapes for atom " << atom << " and element " << el << "
// node "
// << node << " node_global " << node_global
// << " value " << value,
// DBG_ALL);
// auto &m = SubMatrix[el];
// m(atom_index, nd) = value;
// // make the sum shapes for all atoms
// m(m.rows() - 1, nd) += value;
// }
// }
/* --------------------------------------------------------------------------
*/
inline void ShapeMatrix::removeAtom(UInt at, UInt elem) {
// Do not really remove the atom.
// Only substract the contribution to the cumulative part in submatrix.
auto &mat = SubMatrix[elem];
auto &indA = elementToAtoms[elem];
UInt i = 0;
for (i = 0; i < mat.rows() - 1; ++i)
if (indA[i] == at)
break;
LM_ASSERT(i != mat.rows() - 1, "could not find the atom to remove !!");
for (UInt j = 0; j < mat.cols(); ++j) {
mat(mat.rows() - 1, j) -= mat(i, j);
}
indA[i] = UINT_MAX;
}
/* --------------------------------------------------------------------------
*/
inline void ShapeMatrix::buildNodeShape(ContainerArray<Real> &node_shape) {
for (UInt k = 0; k < nb_elements; ++k) {
auto &indN = elementToNodes[k];
auto &mat = SubMatrix[k];
for (UInt j = 0; j < mat.cols(); ++j) {
DUMP("Summing node_shape[" << indN[j] << " (el= " << k << ", j = " << j
<< ")] with " << mat(mat.rows() - 1, j),
DBG_ALL);
LM_ASSERT(mat(mat.rows() - 1, j) >= 0,
"wrong shape fonction in m("
<< mat.rows() - 1 << ", " << j
<< ") = " << &mat(mat.rows() - 1, j) << " = "
<< mat(mat.rows() - 1, j));
node_shape[indN[j]] += mat(mat.rows() - 1, j);
}
}
}
/* --------------------------------------------------------------------------
*/
inline void
ShapeMatrix::buildThermalStuff([[gnu::unused]] std::vector<Real> &nodal_out,
[[gnu::unused]] std::vector<Real> &atom_shape) {
LM_TOIMPLEMENT;
// for (UInt elem = 0; elem < nb_elements; ++elem) {
// UInt *indN = elementToNodes[elem];
// auto &indA = elementToAtoms[elem];
// auto &mat = SubMatrix[elem];
// for (UInt i = 0; i < mat.rows() - 1; ++i) {
// for (UInt j = 0; j < mat.cols(); ++j) {
// nodal_out[indN[j]] += mat(i, j) * atom_shape[indA[i]];
// }
// }
// }
}
/* --------------------------------------------------------------------------
*/
inline void
ShapeMatrix::buildContinuumConstraint(ContainerArray<Real> &A,
ContainerArray<Real> &node_shape,
ContainerArray<Real> &node_weight) {
for (UInt k = 0; k < nb_elements; ++k) {
auto &indA = elementToAtoms[k];
auto &indN = elementToNodes[k];
auto &mat = SubMatrix[k];
for (UInt i = 0; i < mat.rows() - 1; ++i)
if (indA[i] != UINT_MAX)
for (UInt j = 0; j < mat.cols(); ++j) {
DUMP("Building A["
<< indA[i] << "]"
<< " += " << mat(i, j) << "*"
<< " " << node_shape[indN[j]] << " / "
<< node_weight[indN[j]] << std::endl
<< " = "
<< mat(i, j) * node_shape[indN[j]] / node_weight[indN[j]],
DBG_ALL);
A[indA[i]] += mat(i, j) * node_shape[indN[j]] / node_weight[indN[j]];
LM_ASSERT(!std::isinf(A[indA[i]]) && !std::isnan(A[indA[i]]),
"node_weight[" << indN[j] << "] = " << node_weight[indN[j]]
<< "node_shape[" << indN[j]
<< "] = " << node_shape[indN[j]]
<< "fail on atom " << i);
}
}
}
/* --------------------------------------------------------------------------
*/
// template <UInt Dim, typename _Vec_>
// void ShapeMatrix::buildRHS(ContainerArray<Real> &rhs, _Vec_ &v) {
// for (UInt k = 0; k < nb_elements; ++k) {
// auto &mat = SubMatrix[k];
// auto &indGN = elementToNodesGlobalIndex[k];
// auto &indA = elementToAtoms[k];
// // loop over all the atoms of the bloc
// for (UInt i = 0; i < mat.rows() - 1; ++i)
// if (indA[i] != UINT_MAX) {
// VectorView<Dim> _rhs(&rhs(indA[i], 0));
// // summing contributions ovre the nodes
// for (UInt j = 0; j < mat.cols(); ++j) {
// VectorView<Dim> _v(&v(indGN[j], 0));
// DUMP("Before rhs[" << indA[i] << "]=" << _rhs, DBG_ALL);
// _rhs += _v * mat(i, j);
// DUMP("Building " << _rhs << "=rhs[" << indA[i] << "] with "
// << mat(i, j) << "*" << _v << " = " << _v *
// mat(i, j)
// << " (index = " << indGN[j] << ")",
// DBG_ALL);
// }
// }
// }
// }
/* --------------------------------------------------------------------------
*/
template <UInt Dim, typename _Vec_>
void ShapeMatrix::applyContinuumConstraint(_Vec_ &velocity,
ContainerArray<Real> &multL,
ContainerArray<Real> &lambdas) {
// loop over the blocs
for (UInt k = 0; k < nb_elements; ++k) {
auto &mat = SubMatrix[k];
// get indirection to global node indices
// get direct access to array
auto &indGN = elementToNodesGlobalIndex[k];
auto &indN = elementToNodes[k];
auto &indA = elementToAtoms[k];
// loop over the atoms of the bloc
for (UInt i = 0; i < mat.rows() - 1; ++i)
if (indA[i] != UINT_MAX) {
VectorView<Dim> _multL(&multL(indA[i], 0));
// loop over the nodes
for (UInt j = 0; j < mat.cols(); ++j) {
VectorView<Dim> _v(&velocity(indGN[j], 0));
DUMP("adding v[" << indGN[j] << "] with " << -1.0 * mat(i, j) << " * "
<< _multL << " {rhs[" << indA[i] << "]}/ "
<< lambdas[indN[j]] << " = "
<< -1.0 * mat(i, j) * _multL / lambdas[indN[j]],
DBG_ALL);
_v += -1.0 * mat(i, j) * _multL / lambdas[indN[j]];
}
}
}
}
/* --------------------------------------------------------------------------
*/
template <UInt Dim, typename _Vec_>
void ShapeMatrix::applyChangingThermal(_Vec_ &vDisp, _Vec_ &vVel, _Vec_ &vAcc,
_Vec_ &masses, std::vector<Real> &multL,
Real) {
// loop over the blocs
for (UInt k = 0; k < nb_elements; ++k) {
auto &mat = SubMatrix[k];
// get indirection to global node indices
// get direct access to array
auto &indGN = elementToNodesGlobalIndex[k];
auto &indA = elementToAtoms[k];
// loop over the atoms of the block
for (UInt i = 0; i < mat.rows() - 1; ++i)
if (indA[i] != UINT_MAX)
// loop over the nodes
for (UInt j = 0; j < mat.cols(); ++j) {
vDisp.add(indGN[j] * Dim, -1.0 * mat(i, j) * multL[indA[i] * Dim] /
masses[indGN[j] * Dim]);
vDisp.add(indGN[j] * Dim + 1, -1.0 * mat(i, j) *
multL[indA[i] * Dim + 1] /
masses[indGN[j] * Dim + 1]);
vDisp.add(indGN[j] * Dim + 2, -1.0 * mat(i, j) *
multL[indA[i] * Dim + 2] /
masses[indGN[j] * Dim + 2]);
DUMP("Dispx( " << indGN[j] << "," << vDisp[indGN[j] * Dim]
<< ")+=" << mat(i, j) << "*" << multL[indA[i] * Dim]
<< "/" << masses[indGN[j] * Dim],
DBG_ALL);
DUMP("Dispy( " << indGN[j] << "," << vDisp[indGN[j] * Dim + 1]
<< ")+=" << mat(i, j) << "*"
<< multL[indA[i] * Dim + 1] << "/"
<< masses[indGN[j] * Dim + 1],
DBG_ALL);
DUMP("Dispz( " << indGN[j] << "," << vDisp[indGN[j] * Dim + 2]
<< ")+=" << mat(i, j) << "*"
<< multL[indA[i] * Dim + 2] << "/"
<< masses[indGN[j] * Dim + 2],
DBG_ALL);
vVel(indGN[j] * Dim) = 0;
vVel(indGN[j] * Dim + 1) = 0;
vVel(indGN[j] * Dim + 2) = 0;
vAcc(indGN[j] * Dim) = 0;
vAcc(indGN[j] * Dim + 1) = 0;
vAcc(indGN[j] * Dim + 2) = 0;
}
}
}
/* ------------------------------------------------------------------------ */
// template <UInt Dim, typename _Vec_>
// Vector<Dim> ShapeMatrix::interpolate(UInt at, _Vec_ &field, UInt elem) {
// Vector<Dim> res;
// res.setZero();
// LM_ASSERT(elem < nb_elements, "overflow detected");
// auto &indA = elementToAtoms[elem];
// auto &indGN = elementToNodesGlobalIndex[elem];
// auto &mat = SubMatrix[elem];
// DUMP("interpolating for elem " << elem << " atom " << at, DBG_ALL);
// UInt i;
// for (i = 0; i < mat.rows() - 1; ++i)
// if (indA[i] == at)
// break;
// for (UInt j = 0; j < mat.cols(); ++j) {
// for (UInt deg = 0; deg < Dim; ++deg) {
// res[deg] += mat(i, j) * field(indGN[j], deg);
// }
// }
// return res;
// }
/* --------------------------------------------------------------------------
*/
inline Real ShapeMatrix::interpolate(ContainerArray<Real> &atom_data,
ContainerArray<Real> &node_data) {
UInt stride = atom_data.cols();
for (UInt elem = 0; elem < nb_elements; ++elem) {
auto &indN = elementToNodes[elem];
auto &indA = elementToAtoms[elem];
auto &mat = SubMatrix[elem];
for (UInt i = 0; i < mat.rows() - 1; ++i) {
for (UInt j = 0; j < mat.cols(); ++j) {
for (UInt k = 0; k < stride; ++k) {
atom_data(indA[i], k) += mat(i, j) * node_data(indN[j], k);
}
}
}
}
return 0;
}
/* --------------------------------------------------------------------------
*/
inline void ShapeMatrix::getRelatedAtomsToElem(UInt elem,
std::vector<UInt> &atoms) {
auto &indA = elementToAtoms[elem];
auto &mat = SubMatrix[elem];
atoms.clear();
// loop over the atoms
UInt i;
for (i = 0; i < mat.rows() - 1; ++i)
if (indA[i] != UInt(-1))
atoms.push_back(indA[i]);
}
/* --------------------------------------------------------------------------
*/
inline void ShapeMatrix::getRelatedNodesToElem(UInt elem,
std::vector<UInt> &nodes) {
auto &indGN = elementToNodesGlobalIndex[elem];
// auto &mat = SubMatrix[elem];
nodes = indGN;
}
/* --------------------------------------------------------------------------
*/
inline void ShapeMatrix::extractShapeMatrix(Eigen::MatrixXf &A) {
A.setZero();
// loop over elements
for (UInt k = 0; k < nb_elements; ++k) {
auto &indA = elementToAtoms[k];
auto &indN = elementToNodes[k];
auto &mat = SubMatrix[k];
// loop over nodes
for (UInt j = 0; j < mat.cols(); ++j) {
// loop over atoms
for (UInt i = 0; i < mat.rows() - 1; ++i) {
A(indA[i], indN[j]) += mat(i, j);
}
}
}
}
/* --------------------------------------------------------------------------
*/
inline void ShapeMatrix::buildLeastSquareMatrix(Eigen::MatrixXf &A) {
A.setZero();
/* on veut construire les A(h,k) */
/* le produit de deux shape sur la position d'un atome
n'est nul que si les deux shapes font partie
d'un element commun : on a une seul boucle pricipale sur
les elements */
for (UInt k = 0; k < nb_elements; ++k) {
auto &indN = elementToNodes[k];
auto &mat = SubMatrix[k];
/* ici on met la Real boucle sur les noeuds pour calculer
la sous matrice de la grande */
for (UInt j1 = 0; j1 < mat.cols(); ++j1)
for (UInt j2 = 0; j2 < mat.cols(); ++j2) {
/* ici on met la boucle de la somme sur les atomes */
for (UInt i = 0; i < mat.rows() - 1; ++i) {
A(indN[j1], indN[j2]) += mat(i, j1) * mat(i, j2);
DUMP("A[" << indN[j1] << "," << indN[j2]
<< "] = " << A(indN[j1], indN[j2]),
DBG_ALL);
}
}
}
std::ofstream f("least-square-constraint.mat");
f << A;
}
/* --------------------------------------------------------------------------
*/
inline void ShapeMatrix::buildLeastSquareRHS(ContainerArray<Real> &rhs,
ContainerArray<Real> &v) {
UInt stride = rhs.cols();
DUMP("nb_elements: " << nb_elements, DBG_DETAIL);
// loop on elements
for (UInt elem = 0; elem < nb_elements; ++elem) {
DUMP("call build LS RHS. passing elem " << elem, DBG_DETAIL);
auto &indA = elementToAtoms[elem];
auto &indN = elementToNodes[elem];
auto &mat = SubMatrix[elem];
for (UInt i = 0; i < mat.rows() - 1; ++i) // loop on atoms
for (UInt j = 0; j < mat.cols(); ++j) // loop on nodes
for (UInt k = 0; k < stride; ++k) { // loop on dimensions
rhs(indN[j], k) += mat(i, j) * v(indA[i], k);
}
}
}
/* --------------------------------------------------------------------------
*/
inline void ShapeMatrix::buildFormeFaibleCMatrix(Eigen::MatrixXf &C,
std::vector<Real> &lambdas) {
/* le produit de deux shape sur la position d'un atome
n'est nul que si les deux shapes font partie
d'un element commun : on a une seul boucle pricipale sur
les elements */
for (UInt k = 0; k < nb_elements; ++k) {
auto &indN = elementToNodes[k];
auto &indA = elementToAtoms[k];
auto &mat = SubMatrix[k];
/* ici on met la Real boucle sur les noeuds pour calculer
la sous matrice de la grande */
for (UInt j1 = 0; j1 < mat.cols(); ++j1)
for (UInt j2 = 0; j2 < mat.cols(); ++j2) {
/* ici on met la boucle de la somme sur les atomes */
for (UInt i = 0; i < mat.rows() - 1; ++i) {
C(indN[j1], indN[j2]) += mat(i, j1) * mat(i, j2) / lambdas[indA[i]];
DUMP("C[" << indN[j1] << "," << indN[j2]
<< "] = " << C(indN[j1], indN[j2]),
DBG_ALL);
}
}
}
}
/* --------------------------------------------------------------------------
*/
inline void ShapeMatrix::correctFormeFaible(
[[gnu::unused]] Real &v, [[gnu::unused]] Array &rhs,
[[gnu::unused]] std::vector<Real> &lamb, [[gnu::unused]] UInt at,
[[gnu::unused]] UInt elem) {
LM_TOIMPLEMENT;
// UInt *indN = elementToNodes[elem];
// UInt *indA = elementToAtoms[elem];
// auto &mat = SubMatrix[elem];
// /* trouve l'atome concerne */
// UInt i, j;
// for (i = 0; i < mat.rows() - 1; ++i)
// if (indA[i] == at)
// break;
// for (j = 0; j < mat.cols(); ++j) {
// v -= mat(i, j) * rhs(indN[j]) / lamb[at];
// }
}
/* ------------------------------------------------------------------------ */
inline void ShapeMatrix::swapAtomIndirections() {
// UInt **temp = currentElementToAtom;
// currentElementToAtom = newIndirectionAtom;
// newIndirectionAtom = temp;
elementToAtoms.swap();
}
/* ------------------------------------------------------------------------ */
inline void ShapeMatrix::averageOnElements(ContainerArray<Real> &DataA,
ContainerArray<Real> &DataM) {
DataM.setZero();
// loop over elements
for (UInt elem = 0; elem < nb_elements; ++elem) {
auto &indA = elementToAtoms[elem];
auto &mat = SubMatrix[elem];
DataM.row(elem).setZero();
// loop over atoms
for (UInt i = 0; i < mat.rows() - 1; ++i) {
DataM.row(elem) += DataA.row(indA[i]);
}
DataM.row(elem) /= (mat.rows() - 1);
}
}
/* --------------------------------------------------------------------------
*/
inline void ShapeMatrix::alterMatrixIndexesForPBC(
std::vector<std::pair<UInt, UInt>> &pbc_pairs) {
for (UInt k = 0; k < nb_elements; ++k) {
auto &indN = elementToNodes[k];
auto &mat = SubMatrix[k];
for (UInt j = 0; j < mat.cols(); ++j) {
for (UInt p = 0; p < pbc_pairs.size(); p++) {
UInt i1 = pbc_pairs[p].first;
UInt i2 = pbc_pairs[p].second;
if (indN[j] == i1)
indN[j] = i2;
}
}
}
}
/* --------------------------------------------------------------------------
*/
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_SHAPE_MATRIX_HH__ */

Event Timeline

c4science · Help