Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
amdis
amdis-core
Commits
6efe9297
Commit
6efe9297
authored
Mar 20, 2019
by
Praetorius, Simon
Browse files
renamed simplify to as_scalar and add as_vector
parent
0215302f
Changes
2
Hide whitespace changes
Inline
Side-by-side
src/amdis/common/FieldMatVec.hpp
View file @
6efe9297
...
...
@@ -82,51 +82,70 @@ namespace Dune
/// Convert pseudo-scalar to real scalar type
/// @{
template
<
class
T
>
decltype
(
auto
)
simplify
(
T
&&
t
)
decltype
(
auto
)
as_scalar
(
T
&&
t
)
{
return
FWD
(
t
);
}
template
<
class
T
>
T
simplify
(
FieldVector
<
T
,
1
>
const
&
t
)
T
as_scalar
(
FieldVector
<
T
,
1
>
const
&
t
)
{
return
t
[
0
];
}
template
<
class
T
>
T
simplify
(
FieldMatrix
<
T
,
1
,
1
>
const
&
t
)
T
as_scalar
(
FieldMatrix
<
T
,
1
,
1
>
const
&
t
)
{
return
t
[
0
][
0
];
}
template
<
class
T
>
T
simplify
(
DiagonalMatrix
<
T
,
1
>
const
&
t
)
T
as_scalar
(
DiagonalMatrix
<
T
,
1
>
const
&
t
)
{
return
t
.
diagonal
(
0
);
}
/// @}
/// Convert pseudo-vector to real vector type
/// @{
template
<
class
T
>
decltype
(
auto
)
as_vector
(
T
&&
t
)
{
return
FWD
(
t
);
}
template
<
class
T
,
int
N
>
FieldVector
<
T
,
N
>
const
&
as_vector
(
FieldMatrix
<
T
,
1
,
N
>
const
&
t
)
{
return
t
[
0
];
}
template
<
class
T
,
int
N
>
FieldVector
<
T
,
N
>&
as_vector
(
FieldMatrix
<
T
,
1
,
N
>&
t
)
{
return
t
[
0
];
}
/// @}
// returns -a
template
<
class
A
>
auto
negate
(
A
const
&
a
);
// returns a+b
template
<
class
A
,
class
B
>
auto
plus
(
A
a
,
B
const
&
b
)
{
return
a
+=
b
;
}
auto
plus
(
A
const
&
a
,
B
const
&
b
);
// returns a-b
template
<
class
A
,
class
B
>
auto
minus
(
A
a
,
B
const
&
b
)
{
return
a
-=
b
;
}
auto
minus
(
A
const
&
a
,
B
const
&
b
);
// returns a*b
template
<
class
A
,
class
B
,
std
::
enable_if_t
<
IsNumber
<
A
>
::
value
||
IsNumber
<
B
>::
value
,
int
>
=
0
>
std
::
enable_if_t
<
IsNumber
<
A
>
::
value
&&
IsNumber
<
B
>::
value
,
int
>
=
0
>
auto
multiplies
(
A
const
&
a
,
B
const
&
b
);
template
<
class
A
,
class
B
,
std
::
enable_if_t
<
IsNumber
<
A
>
::
value
!=
IsNumber
<
B
>::
value
,
int
>
=
0
>
auto
multiplies
(
A
const
&
a
,
B
const
&
b
);
template
<
class
T
,
int
N
,
class
S
>
...
...
@@ -158,35 +177,35 @@ namespace Dune
std
::
enable_if_t
<
MatVec
::
IsMatVec
<
A
>
::
value
,
int
>
=
0
>
auto
operator
-
(
A
const
&
a
)
{
return
MatVec
::
negate
(
MatVec
::
simplify
(
a
));
return
MatVec
::
negate
(
MatVec
::
as_scalar
(
a
));
}
template
<
class
A
,
class
B
,
std
::
enable_if_t
<
MatVec
::
IsMatVec
<
A
>
::
value
||
MatVec
::
IsMatVec
<
B
>::
value
,
int
>
=
0
>
auto
operator
+
(
A
const
&
a
,
B
const
&
b
)
{
return
MatVec
::
plus
(
MatVec
::
simplify
(
a
),
MatVec
::
simplify
(
b
));
return
MatVec
::
plus
(
MatVec
::
as_scalar
(
a
),
MatVec
::
as_scalar
(
b
));
}
template
<
class
A
,
class
B
,
std
::
enable_if_t
<
MatVec
::
IsMatVec
<
A
>
::
value
||
MatVec
::
IsMatVec
<
B
>::
value
,
int
>
=
0
>
auto
operator
-
(
A
const
&
a
,
B
const
&
b
)
{
return
MatVec
::
minus
(
MatVec
::
simplify
(
a
),
MatVec
::
simplify
(
b
));
return
MatVec
::
minus
(
MatVec
::
as_scalar
(
a
),
MatVec
::
as_scalar
(
b
));
}
template
<
class
A
,
class
B
,
std
::
enable_if_t
<
MatVec
::
IsMatVec
<
A
>
::
value
||
MatVec
::
IsMatVec
<
B
>::
value
,
int
>
=
0
>
auto
operator
*
(
A
const
&
a
,
B
const
&
b
)
{
return
MatVec
::
multiplies
(
MatVec
::
simplify
(
a
),
MatVec
::
simplify
(
b
));
return
MatVec
::
multiplies
(
MatVec
::
as_scalar
(
a
),
MatVec
::
as_scalar
(
b
));
}
template
<
class
A
,
class
B
,
std
::
enable_if_t
<
MatVec
::
IsMatVec
<
A
>
::
value
||
MatVec
::
IsMatVec
<
B
>::
value
,
int
>
=
0
>
auto
operator
/
(
A
const
&
a
,
B
const
&
b
)
{
return
MatVec
::
divides
(
MatVec
::
simplify
(
a
),
MatVec
::
simplify
(
b
));
return
MatVec
::
divides
(
MatVec
::
as_scalar
(
a
),
MatVec
::
as_scalar
(
b
));
}
// ----------------------------------------------------------------------------
...
...
@@ -346,19 +365,28 @@ namespace Dune
template
<
class
T
,
int
M
,
int
N
>
FieldMatrix
<
T
,
N
,
M
>
trans
(
FieldMatrix
<
T
,
M
,
N
>
const
&
A
);
template
<
class
T
,
int
N
>
DiagonalMatrix
<
T
,
N
>
const
&
trans
(
DiagonalMatrix
<
T
,
N
>
const
&
A
)
{
return
A
;
}
// -----------------------------------------------------------------------------
template
<
class
T
,
int
M
,
int
N
,
int
L
>
FieldMatrix
<
T
,
M
,
N
>
multiplies_AtB
(
FieldMatrix
<
T
,
L
,
M
>
const
&
A
,
FieldMatrix
<
T
,
N
,
L
>
const
&
B
);
template
<
class
T
1
,
class
T2
,
int
M
,
int
N
,
int
L
>
FieldMatrix
<
std
::
common_type_t
<
T1
,
T2
>
,
M
,
N
>
multiplies_AtB
(
FieldMatrix
<
T
1
,
L
,
M
>
const
&
A
,
FieldMatrix
<
T
2
,
N
,
L
>
const
&
B
);
template
<
class
T
,
int
M
,
int
N
,
int
L
>
FieldMatrix
<
T
,
M
,
N
>
multiplies_ABt
(
FieldMatrix
<
T
,
M
,
L
>
const
&
A
,
FieldMatrix
<
T
,
N
,
L
>
const
&
B
);
template
<
class
T
1
,
class
T2
,
int
M
,
int
N
,
int
L
>
FieldMatrix
<
std
::
common_type_t
<
T1
,
T2
>
,
M
,
N
>
multiplies_ABt
(
FieldMatrix
<
T
1
,
M
,
L
>
const
&
A
,
FieldMatrix
<
T
2
,
N
,
L
>
const
&
B
);
template
<
class
T
,
int
M
,
int
N
,
int
L
>
FieldMatrix
<
T
,
M
,
N
>&
multiplies_ABt
(
FieldMatrix
<
T
,
M
,
L
>
const
&
A
,
FieldMatrix
<
T
,
N
,
L
>
const
&
B
,
FieldMatrix
<
T
,
M
,
N
>&
C
);
template
<
class
T
1
,
class
T2
,
class
T3
,
int
M
,
int
N
,
int
L
>
FieldMatrix
<
T
3
,
M
,
N
>&
multiplies_ABt
(
FieldMatrix
<
T
1
,
M
,
L
>
const
&
A
,
FieldMatrix
<
T
2
,
N
,
L
>
const
&
B
,
FieldMatrix
<
T
3
,
M
,
N
>&
C
);
template
<
class
T
,
int
M
,
int
N
>
FieldMatrix
<
T
,
M
,
N
>&
multiplies_ABt
(
FieldMatrix
<
T
,
M
,
N
>
const
&
A
,
DiagonalMatrix
<
T
,
N
>
const
&
B
,
FieldMatrix
<
T
,
M
,
N
>&
C
);
template
<
class
T1
,
class
T2
,
class
T3
,
int
M
,
int
N
>
FieldMatrix
<
T3
,
M
,
N
>&
multiplies_ABt
(
FieldMatrix
<
T1
,
M
,
N
>
const
&
A
,
DiagonalMatrix
<
T2
,
N
>
const
&
B
,
FieldMatrix
<
T3
,
M
,
N
>&
C
);
template
<
class
T1
,
class
T2
,
class
T3
,
int
N
>
FieldVector
<
T3
,
N
>&
multiplies_ABt
(
FieldMatrix
<
T1
,
1
,
N
>
const
&
A
,
DiagonalMatrix
<
T2
,
N
>
const
&
B
,
FieldVector
<
T3
,
N
>&
C
);
// -----------------------------------------------------------------------------
...
...
src/amdis/common/FieldMatVec.inc.hpp
View file @
6efe9297
...
...
@@ -3,6 +3,8 @@
#include <algorithm>
#include <limits>
#include <dune/functions/functionspacebases/flatvectorview.hh>
#include <amdis/common/Math.hpp>
#include <amdis/operations/Arithmetic.hpp>
#include <amdis/operations/MaxMin.hpp>
...
...
@@ -19,8 +21,47 @@ namespace MatVec {
return
multiplies
(
a
,
-
1
);
}
// returns a+b
template
<
class
A
,
class
B
>
auto
plus
(
A
const
&
a
,
B
const
&
b
)
{
using
T
=
std
::
common_type_t
<
typename
FieldTraits
<
A
>::
field_type
,
typename
FieldTraits
<
B
>::
field_type
>
;
typename
MakeMatVec
<
A
,
T
>::
type
c
{
a
};
auto
b_
=
Dune
::
Functions
::
flatVectorView
(
b
);
auto
c_
=
Dune
::
Functions
::
flatVectorView
(
c
);
assert
(
int
(
b_
.
size
())
==
int
(
c_
.
size
()));
for
(
int
i
=
0
;
i
<
int
(
b_
.
size
());
++
i
)
c_
[
i
]
+=
b_
[
i
];
return
c
;
}
// returns a-b
template
<
class
A
,
class
B
>
auto
minus
(
A
const
&
a
,
B
const
&
b
)
{
using
T
=
std
::
common_type_t
<
typename
FieldTraits
<
A
>::
field_type
,
typename
FieldTraits
<
B
>::
field_type
>
;
typename
MakeMatVec
<
A
,
T
>::
type
c
{
a
};
auto
b_
=
Dune
::
Functions
::
flatVectorView
(
b
);
auto
c_
=
Dune
::
Functions
::
flatVectorView
(
c
);
assert
(
int
(
b_
.
size
())
==
int
(
c_
.
size
()));
for
(
int
i
=
0
;
i
<
int
(
b_
.
size
());
++
i
)
c_
[
i
]
-=
b_
[
i
];
return
c
;
}
template
<
class
A
,
class
B
,
std
::
enable_if_t
<
IsNumber
<
A
>
::
value
||
IsNumber
<
B
>::
value
,
int
>>
std
::
enable_if_t
<
IsNumber
<
A
>
::
value
&&
IsNumber
<
B
>::
value
,
int
>>
auto
multiplies
(
A
const
&
a
,
B
const
&
b
)
{
return
a
*
b
;
}
template
<
class
A
,
class
B
,
std
::
enable_if_t
<
IsNumber
<
A
>
::
value
!=
IsNumber
<
B
>::
value
,
int
>>
auto
multiplies
(
A
const
&
a
,
B
const
&
b
)
{
using
T
=
std
::
common_type_t
<
typename
FieldTraits
<
A
>::
field_type
,
typename
FieldTraits
<
B
>::
field_type
>
;
...
...
@@ -158,14 +199,14 @@ T sum(FieldMatrix<T, 1, N> const& x)
template
<
class
T
,
int
N
>
auto
unary_dot
(
FieldVector
<
T
,
N
>
const
&
x
)
{
auto
op
=
[](
auto
const
&
a
,
auto
const
&
b
)
{
return
a
+
AMDiS
::
Math
::
sqr
(
std
::
abs
(
b
));
};
auto
op
=
[](
auto
const
&
a
,
auto
const
&
b
)
{
using
std
::
abs
;
return
a
+
AMDiS
::
Math
::
sqr
(
abs
(
b
));
};
return
Impl
::
accumulate
(
x
,
T
(
0
),
op
);
}
template
<
class
T
,
int
N
>
auto
unary_dot
(
FieldMatrix
<
T
,
1
,
N
>
const
&
x
)
{
auto
op
=
[](
auto
const
&
a
,
auto
const
&
b
)
{
return
a
+
AMDiS
::
Math
::
sqr
(
std
::
abs
(
b
));
};
auto
op
=
[](
auto
const
&
a
,
auto
const
&
b
)
{
using
std
::
abs
;
return
a
+
AMDiS
::
Math
::
sqr
(
abs
(
b
));
};
return
Impl
::
accumulate
(
x
,
T
(
0
),
op
);
}
...
...
@@ -229,14 +270,14 @@ auto abs_min(FieldMatrix<T, 1, N> const& x)
template
<
class
T
,
int
N
>
auto
one_norm
(
FieldVector
<
T
,
N
>
const
&
x
)
{
auto
op
=
[](
auto
const
&
a
,
auto
const
&
b
)
{
return
a
+
std
::
abs
(
b
);
};
auto
op
=
[](
auto
const
&
a
,
auto
const
&
b
)
{
using
std
::
abs
;
return
a
+
abs
(
b
);
};
return
Impl
::
accumulate
(
x
,
T
(
0
),
op
);
}
template
<
class
T
,
int
N
>
auto
one_norm
(
FieldMatrix
<
T
,
1
,
N
>
const
&
x
)
{
auto
op
=
[](
auto
const
&
a
,
auto
const
&
b
)
{
return
a
+
std
::
abs
(
b
);
};
auto
op
=
[](
auto
const
&
a
,
auto
const
&
b
)
{
using
std
::
abs
;
return
a
+
abs
(
b
);
};
return
Impl
::
accumulate
(
x
,
T
(
0
),
op
);
}
...
...
@@ -261,14 +302,14 @@ auto two_norm(FieldMatrix<T, 1, N> const& x)
template
<
int
p
,
class
T
,
int
N
>
auto
p_norm
(
FieldVector
<
T
,
N
>
const
&
x
)
{
auto
op
=
[](
auto
const
&
a
,
auto
const
&
b
)
{
return
a
+
AMDiS
::
Math
::
pow
<
p
>
(
std
::
abs
(
b
));
};
auto
op
=
[](
auto
const
&
a
,
auto
const
&
b
)
{
using
std
::
abs
;
return
a
+
AMDiS
::
Math
::
pow
<
p
>
(
abs
(
b
));
};
return
std
::
pow
(
Impl
::
accumulate
(
x
,
T
(
0
),
op
),
1.0
/
p
);
}
template
<
int
p
,
class
T
,
int
N
>
auto
p_norm
(
FieldMatrix
<
T
,
1
,
N
>
const
&
x
)
{
auto
op
=
[](
auto
const
&
a
,
auto
const
&
b
)
{
return
a
+
AMDiS
::
Math
::
pow
<
p
>
(
std
::
abs
(
b
));
};
auto
op
=
[](
auto
const
&
a
,
auto
const
&
b
)
{
using
std
::
abs
;
return
a
+
AMDiS
::
Math
::
pow
<
p
>
(
abs
(
b
));
};
return
std
::
pow
(
Impl
::
accumulate
(
x
,
T
(
0
),
op
),
1.0
/
p
);
}
...
...
@@ -293,10 +334,11 @@ auto infty_norm(FieldMatrix<T, 1, N> const& x)
template
<
class
T
,
int
N
>
T
distance
(
FieldVector
<
T
,
N
>
const
&
lhs
,
FieldVector
<
T
,
N
>
const
&
rhs
)
{
using
std
::
sqrt
;
T
result
=
0
;
for
(
int
i
=
0
;
i
<
N
;
++
i
)
result
+=
AMDiS
::
Math
::
sqr
(
lhs
[
i
]
-
rhs
[
i
]);
return
std
::
sqrt
(
result
);
return
sqrt
(
result
);
}
// ----------------------------------------------------------------------------
...
...
@@ -397,10 +439,10 @@ FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A)
}
template
<
class
T
,
int
M
,
int
N
,
int
L
>
FieldMatrix
<
T
,
M
,
N
>
multiplies_AtB
(
FieldMatrix
<
T
,
L
,
M
>
const
&
A
,
FieldMatrix
<
T
,
N
,
L
>
const
&
B
)
template
<
class
T
1
,
class
T2
,
int
M
,
int
N
,
int
L
>
FieldMatrix
<
std
::
common_type_t
<
T1
,
T2
>
,
M
,
N
>
multiplies_AtB
(
FieldMatrix
<
T
1
,
L
,
M
>
const
&
A
,
FieldMatrix
<
T
2
,
N
,
L
>
const
&
B
)
{
FieldMatrix
<
T
,
M
,
N
>
C
;
FieldMatrix
<
std
::
common_type_t
<
T1
,
T2
>
,
M
,
N
>
C
;
for
(
int
m
=
0
;
m
<
M
;
++
m
)
{
for
(
int
n
=
0
;
n
<
N
;
++
n
)
{
...
...
@@ -412,15 +454,15 @@ FieldMatrix<T,M,N> multiplies_AtB(FieldMatrix<T, L, M> const& A, FieldMatrix<T,
return
C
;
}
template
<
class
T
,
int
M
,
int
N
,
int
L
>
FieldMatrix
<
T
,
M
,
N
>
multiplies_ABt
(
FieldMatrix
<
T
,
M
,
L
>
const
&
A
,
FieldMatrix
<
T
,
N
,
L
>
const
&
B
)
template
<
class
T
1
,
class
T2
,
int
M
,
int
N
,
int
L
>
FieldMatrix
<
std
::
common_type_t
<
T1
,
T2
>
,
M
,
N
>
multiplies_ABt
(
FieldMatrix
<
T
1
,
M
,
L
>
const
&
A
,
FieldMatrix
<
T
2
,
N
,
L
>
const
&
B
)
{
FieldMatrix
<
T
,
M
,
N
>
C
;
FieldMatrix
<
std
::
common_type_t
<
T1
,
T2
>
,
M
,
N
>
C
;
return
multiplies_ABt
(
A
,
B
,
C
);
}
template
<
class
T
,
int
M
,
int
N
,
int
L
>
FieldMatrix
<
T
,
M
,
N
>&
multiplies_ABt
(
FieldMatrix
<
T
,
M
,
L
>
const
&
A
,
FieldMatrix
<
T
,
N
,
L
>
const
&
B
,
FieldMatrix
<
T
,
M
,
N
>&
C
)
template
<
class
T
1
,
class
T2
,
class
T3
,
int
M
,
int
N
,
int
L
>
FieldMatrix
<
T
3
,
M
,
N
>&
multiplies_ABt
(
FieldMatrix
<
T
1
,
M
,
L
>
const
&
A
,
FieldMatrix
<
T
2
,
N
,
L
>
const
&
B
,
FieldMatrix
<
T
3
,
M
,
N
>&
C
)
{
for
(
int
m
=
0
;
m
<
M
;
++
m
)
{
for
(
int
n
=
0
;
n
<
N
;
++
n
)
{
...
...
@@ -432,8 +474,8 @@ FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, L> const& A, FieldMatrix<T
return
C
;
}
template
<
class
T
,
int
M
,
int
N
>
FieldMatrix
<
T
,
M
,
N
>&
multiplies_ABt
(
FieldMatrix
<
T
,
M
,
N
>
const
&
A
,
DiagonalMatrix
<
T
,
N
>
const
&
B
,
FieldMatrix
<
T
,
M
,
N
>&
C
)
template
<
class
T
1
,
class
T2
,
class
T3
,
int
M
,
int
N
>
FieldMatrix
<
T
3
,
M
,
N
>&
multiplies_ABt
(
FieldMatrix
<
T
1
,
M
,
N
>
const
&
A
,
DiagonalMatrix
<
T
2
,
N
>
const
&
B
,
FieldMatrix
<
T
3
,
M
,
N
>&
C
)
{
for
(
int
m
=
0
;
m
<
M
;
++
m
)
{
for
(
int
n
=
0
;
n
<
N
;
++
n
)
{
...
...
@@ -443,6 +485,15 @@ FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, N> const& A, DiagonalMatri
return
C
;
}
template
<
class
T1
,
class
T2
,
class
T3
,
int
N
>
FieldVector
<
T3
,
N
>&
multiplies_ABt
(
FieldMatrix
<
T1
,
1
,
N
>
const
&
A
,
DiagonalMatrix
<
T2
,
N
>
const
&
B
,
FieldVector
<
T3
,
N
>&
C
)
{
for
(
int
n
=
0
;
n
<
N
;
++
n
)
{
C
[
n
]
=
A
[
0
][
n
]
*
B
.
diagonal
(
n
);
}
return
C
;
}
template
<
class
T
,
int
N
>
T
const
&
at
(
FieldMatrix
<
T
,
N
,
1
>
const
&
vec
,
std
::
size_t
i
)
{
...
...
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment