Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F82603303
shape_matrix.hh
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Sep 12, 09:12
Size
27 KB
Mime Type
text/x-c++
Expires
Sat, Sep 14, 09:12 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
20721001
Attached To
rLIBMULTISCALE LibMultiScale
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__ */
Event Timeline
Log In to Comment