Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F87420603
system_base.cc
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
Sat, Oct 12, 14:23
Size
6 KB
Mime Type
text/x-c
Expires
Mon, Oct 14, 14:23 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
21587690
Attached To
rMSPPROTO µSpectre prototype implementation
system_base.cc
View Options
/**
* file system_base.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 09 May 2017
*
* @brief Implementation for system base class
*
* @section LICENCE
*
* Copyright (C) 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include <stdexcept>
#include <algorithm>
#include <utility>
#include "system/system_base.hh"
#include "iostream"
namespace
muSpectre
{
namespace
Projection
{
//----------------------------------------------------------------------------//
//! computes fft frequecies. see function fftfreq in numpy's helper.py
vector
fft_freqs
(
Uint
nb_samples
)
{
vector
retval
(
1
,
nb_samples
);
int
N
=
(
nb_samples
-
1
)
/
2
+
1
;
// needs to be signed int for neg freqs
for
(
int
i
=
0
;
i
<
N
;
++
i
)
{
retval
(
0
,
i
)
=
i
;
}
for
(
int
i
=
N
;
i
<
static_cast
<
int
>
(
nb_samples
);
++
i
)
{
retval
(
0
,
i
)
=
-
(
static_cast
<
int
>
(
nb_samples
)
/
2
)
+
i
-
N
;
}
return
retval
;
}
//----------------------------------------------------------------------------//
vector
fft_freqs
(
Uint
nb_samples
,
Real
length
)
{
return
fft_freqs
(
nb_samples
)
/
length
;
}
//----------------------------------------------------------------------------//
template
<
Dim_t
Dim
>
FreqStruc
<
Dim
>::
FreqStruc
(
const
std
::
array
<
Real
,
Dim
>&
sizes
,
const
Index_t
<
Dim
>&
nb_pixels
)
:
nb_pixels
(
nb_pixels
),
sizes
(
sizes
)
{
for
(
size_t
i
=
0
;
i
<
Dim
;
++
i
)
{
this
->
xi
[
i
]
=
Projection
::
fft_freqs
(
this
->
nb_pixels
[
i
]);
}
}
//----------------------------------------------------------------------------//
template
<
Dim_t
Dim
>
inline
typename
FreqStruc
<
Dim
>::
iterator
FreqStruc
<
Dim
>::
begin
()
{
return
iterator
(
*
this
,
{
0
});
}
//----------------------------------------------------------------------------//
template
<
Dim_t
Dim
>
inline
typename
FreqStruc
<
Dim
>::
iterator
FreqStruc
<
Dim
>::
end
()
{
return
iterator
(
*
this
,
this
->
nb_pixels
);
}
//----------------------------------------------------------------------------//
template
<
Dim_t
Dim
>
FreqStruc
<
Dim
>::
iterator
::
iterator
(
const
FreqStruc
<
Dim
>
&
parent
,
Index_t
<
Dim
>
index
)
:
parent
(
parent
),
index
(
std
::
accumulate
(
index
.
begin
(),
index
.
end
(),
1
,
std
::
multiplies
<
size_t
>
())){
}
//----------------------------------------------------------------------------//
template
<
Dim_t
Dim
>
inline
Index_t
<
Dim
>
FreqStruc
<
Dim
>::
iterator
::
get_index
()
{
// implements row-major choice. might have to be templated at some point
static_assert
(((
Dim
==
twoD
)
||
(
Dim
==
threeD
)),
"only 2d or 3d"
);
auto
&
stride
=
this
->
parent
.
nb_pixels
[
1
];
return
Index_t
<
Dim
>
{
static_cast
<
Dim_t
>
(
this
->
index
/
stride
),
static_cast
<
Dim_t
>
(
this
->
index
%
stride
)};
}
//----------------------------------------------------------------------------//
template
<>
inline
Index_t
<
threeD
>
FreqStruc
<
threeD
>::
iterator
::
get_index
()
{
// implements row-major choice. might have to be templated at some point
auto
&&
kstride
=
this
->
parent
.
nb_pixels
[
2
];
Dim_t
&&
k
=
this
->
index
%
kstride
;
auto
&&
jstride
=
this
->
parent
.
nb_pixels
[
1
];
Dim_t
&&
j
=
(
this
->
index
/
kstride
)
%
jstride
;
Dim_t
&&
i
=
(
this
->
index
/
kstride
)
/
jstride
;
return
Index_t
<
threeD
>
{
i
,
j
,
k
};
}
//----------------------------------------------------------------------------//
template
<
Dim_t
Dim
>
inline
Eigen
::
Matrix
<
Real
,
Dim
,
1
>
FreqStruc
<
Dim
>::
iterator
::
operator
*
()
{
Eigen
::
Matrix
<
Real
,
Dim
,
1
>
retval
;
auto
&&
indices
=
this
->
get_index
();
if
(
Dim
==
twoD
)
{
retval
<<
parent
.
xi
[
0
](
indices
[
0
]),
parent
.
xi
[
1
](
indices
[
1
]);
}
else
if
(
Dim
==
threeD
)
{
retval
<<
parent
.
xi
[
0
](
indices
[
0
]),
parent
.
xi
[
1
](
indices
[
1
]),
parent
.
xi
[
2
](
indices
[
2
]);
}
return
retval
;
}
//----------------------------------------------------------------------------//
template
<
Dim_t
Dim
>
inline
typename
FreqStruc
<
Dim
>::
iterator
&
FreqStruc
<
Dim
>::
iterator
::
operator
++
()
{
this
->
index
++
;
return
*
this
;
}
}
// Projection
//----------------------------------------------------------------------------//
template
<
Dim_t
DimS
,
Dim_t
DimM
>
SystemBase
<
DimS
,
DimM
>::
SystemBase
(
std
::
array
<
Real
,
DimS
>
sizes
,
Index
nb_pixels
,
Formulation
form
)
:
sizes
{
sizes
},
nb_pixels
{
nb_pixels
},
nb_pixel
(
std
::
accumulate
(
nb_pixels
.
begin
(),
nb_pixels
.
end
(),
1
,
std
::
multiplies
<
size_t
>
())),
form
{
form
}
{
for
(
const
auto
&
size:
this
->
nb_pixels
)
{
if
(
size
%
2
!=
1
)
{
throw
std
::
runtime_error
(
"Sorry, for now can only handle odd grid sizes"
);
}
}
if
(
form
!=
Formulation
::
finite_strain
)
{
throw
std
::
runtime_error
(
"Sorry, only finite strain for the time being"
);
}
//resize and allocate
this
->
Ghats
=
T4Vector
(
DimM
,
DimM
,
DimM
,
DimM
,
this
->
size
());
this
->
Ghats
.
setZero
();
this
->
build_ghats
();
}
//----------------------------------------------------------------------------//
template
<
Dim_t
DimS
,
Dim_t
DimM
>
void
SystemBase
<
DimS
,
DimM
>::
build_ghats
()
{
//! row-major storage for pixels
Projection
::
FreqStruc
<
DimS
>
freqs
(
this
->
sizes
,
this
->
nb_pixels
);
auto
it
=
freqs
.
begin
();
//drop the first frequency, as it is zero
++
it
;
for
(
size_t
pix_id
=
1
;
pix_id
<
this
->
nb_pixel
;
++
pix_id
,
++
it
)
{
auto
xi
=
*
it
;
xi
/=
xi
.
norm
();
using
T4shape
=
Eigen
::
Sizes
<
DimS
,
DimS
,
DimS
,
DimS
>
;
using
T4
=
Eigen
::
TensorFixedSize
<
Real
,
T4shape
>
;
Eigen
::
TensorMap
<
T4
>
Ghat_kk
(
&
this
->
Ghats
(
0
,
0
,
0
,
0
,
pix_id
),
DimS
,
DimS
,
DimS
,
DimS
);
for
(
Dim_t
im
=
0
;
im
<
DimS
;
++
im
)
{
for
(
Dim_t
j
=
0
;
j
<
DimS
;
++
j
)
{
for
(
Dim_t
l
=
0
;
l
<
DimS
;
++
l
)
{
Ghat_kk
(
im
,
j
,
l
,
im
)
=
xi
(
j
)
*
xi
(
l
);
}
}
}
}
}
template
class
SystemBase
<
2
,
2
>
;
template
class
SystemBase
<
2
,
3
>
;
template
class
SystemBase
<
3
,
3
>
;
}
// muSpectre
Event Timeline
Log In to Comment