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
fa24f5f8
Commit
fa24f5f8
authored
Jan 02, 2018
by
Praetorius, Simon
Browse files
FileWriter cleaned up
parent
ce0fece2
Changes
14
Hide whitespace changes
Inline
Side-by-side
doc/Mainpage.md
View file @
fa24f5f8
...
...
@@ -29,7 +29,7 @@ The corresponding weak form of the equation reads:
with
\f
$ u
\i
n V_g
\f
$.
Thus, we need to define a grid (discretization of
\f
$
\O
mega
\f
$) and a finite
element space (discretization of
\f
$ V_0
\f
$ and
\f
$ V_g
\f
$). This
needs to
be
element space (discretization of
\f
$ V_0
\f
$ and
\f
$ V_g
\f
$). This
cab
be
provided as
`Traits`
parameter in the ProblemStat class:
~~~~~~~~~~~~~~~
{.cpp}
...
...
@@ -37,7 +37,7 @@ using Grid = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>;
using Traits = LagrangeBasis<Grid::LeafGridView, 1>;
~~~~~~~~~~~~~~~
This defines and
`AlbertaGrid`
as grid and
a Lagrange Finite-Element space with
where
`AlbertaGrid`
defines a grid and
`Traits`
a Lagrange Finite-Element space with
local polynomial degree 1 on the elements.
All AMDiS programs start with initialization of the library, using
`AMDiS::init`
...
...
@@ -62,7 +62,7 @@ prob.initialize(INIT_ALL);
Operators specify the (bi-)linear-form and the coefficient function in the term,
see
\r
ef operators for a list of possible types. The bilinear-form in the Poisson
equation consists of second-order term with coefficit
n
= 1 and the linear-form
equation consists of
a
second-order term with coeffici
en
t = 1 and the linear-form
includes the function
\f
$ f(x)=-1
\f
$:
~~~~~~~~~~~~~~~
{.cpp}
...
...
@@ -74,8 +74,8 @@ prob.addVectorOperator(opF, 0);
~~~~~~~~~~~~~~~
Boundary conditions, in the example above a Dirichlet condition, is specified by
defining a predicate for the boundary
\f
$
\p
artial
\O
mega
\f
$ and the
values on the
boundary
\f
$ g(x) = 0
\f
$:
defining a predicate for the boundary
\f
$
\G
amma
\s
ubset
\p
artial
\O
mega
\f
$ and the
values on the
boundary
\f
$ g(x) = 0
\f
$:
~~~~~~~~~~~~~~~
{.cpp}
auto predicate = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8; };
...
...
@@ -93,7 +93,7 @@ adapt.adapt(); // assemble and solve
Finally, finish the AMDiS program with
`AMDiS::finish()`
.
In
complete program then reads:
The
complete program then reads:
~~~~~~~~~~~~~~~
{.cpp}
#include <dune/amdis/AMDiS.hpp>
#include <dune/amdis/AdaptInfo.hpp>
...
...
dune/amdis/FileWriter.hpp
View file @
fa24f5f8
...
...
@@ -10,20 +10,48 @@
#include <dune/typetree/childextraction.hh>
#include <dune/amdis/Initfile.hpp>
#include <dune/amdis/common/Size.hpp>
#include <dune/amdis/common/ValueCategory.hpp>
#include <dune/amdis/gridfunctions/DOFVectorView.hpp>
#include <dune/amdis/io/FileWriterInterface.hpp>
#include <dune/amdis/utility/Filesystem.hpp>
namespace
AMDiS
{
namespace
Impl
{
template
<
class
Tag
>
struct
VTKFieldTypeImpl
;
template
<
>
struct
VTKFieldTypeImpl
<
tag
::
scalar
>
{
static
const
Dune
::
VTK
::
FieldInfo
::
Type
value
=
Dune
::
VTK
::
FieldInfo
::
Type
::
scalar
;
};
template
<
>
struct
VTKFieldTypeImpl
<
tag
::
vector
>
{
static
const
Dune
::
VTK
::
FieldInfo
::
Type
value
=
Dune
::
VTK
::
FieldInfo
::
Type
::
vector
;
};
template
<
>
struct
VTKFieldTypeImpl
<
tag
::
matrix
>
{
static
const
Dune
::
VTK
::
FieldInfo
::
Type
value
=
Dune
::
VTK
::
FieldInfo
::
Type
::
tensor
;
};
}
// end namespace Impl
template
<
class
Range
>
constexpr
Dune
::
VTK
::
FieldInfo
::
Type
VTKFieldType
=
Impl
::
VTKFieldTypeImpl
<
ValueCategory_t
<
Range
>>::
value
;
template
<
class
Range
>
constexpr
std
::
size_t
VTKFieldSize
=
Size
<
Range
>
;
template
<
class
Traits
,
class
TreePath
>
class
FileWriter
:
public
FileWriterInterface
{
private:
// typedefs and static constants
using
GlobalBasis
=
typename
Traits
::
GlobalBasis
;
using
GridView
=
typename
GlobalBasis
::
GridView
;
using
Vector
=
DOFVectorConstView
<
GlobalBasis
,
TreePath
>
;
using
Range
=
typename
Vector
::
Range
;
/// Dimension of the mesh
static
constexpr
int
dim
=
GridView
::
dimension
;
...
...
@@ -31,117 +59,57 @@ namespace AMDiS
/// Dimension of the world
static
constexpr
int
dow
=
GridView
::
dimensionworld
;
using
Vector
=
DOFVector
<
GlobalBasis
>
;
public:
/// Constructor.
FileWriter
(
std
::
string
baseName
,
std
::
shared_ptr
<
GlobalBasis
>
const
&
basis
,
TreePath
const
&
tp
,
std
::
shared_ptr
<
Vector
>
const
&
vector
)
Vector
const
&
dofvector
)
:
FileWriterInterface
(
baseName
)
,
basis_
(
basis
)
,
treePath_
(
tp
)
,
vector_
(
vector
)
,
dofvector_
(
dofvector
)
{
//int subSampling = Parameters::get<int>(baseName + "->subsampling").value_or(0);
int
m
=
Parameters
::
get
<
int
>
(
baseName
+
"->ParaView mode"
).
value_or
(
0
);
mode_
=
(
m
==
0
?
Dune
::
VTK
::
ascii
:
Dune
::
VTK
::
appendedraw
);
//if (subSampling > 0)
// vtkWriter_ = std::make_shared<Dune::SubsamplingVTKWriter<GridView>>(basis_->gridView(), subSampling);
//else
vtkWriter_
=
std
::
make_shared
<
Dune
::
VTKWriter
<
GridView
>>
(
basis_
->
gridView
());
int
subSampling
=
Parameters
::
get
<
int
>
(
baseName
+
"->subsampling"
).
value_or
(
0
);
vtkSeqWriter_
=
std
::
make_shared
<
Dune
::
VTKSequenceWriter
<
GridView
>>
(
vtkWriter_
,
filename_
,
dir_
,
""
);
if
(
subSampling
>
0
)
vtkWriter_
=
std
::
make_shared
<
Dune
::
SubsamplingVTKWriter
<
GridView
>>
(
gridView
(),
subSampling
);
else
vtkWriter_
=
std
::
make_shared
<
Dune
::
VTKWriter
<
GridView
>>
(
gridView
());
mode_
=
Parameters
::
get
<
int
>
(
baseName
+
"->ParaView mode"
).
value_or
(
0
);
vtkSeqWriter_
=
std
::
make_shared
<
Dune
::
VTKSequenceWriter
<
GridView
>>
(
vtkWriter_
,
filename_
,
dir_
,
""
);
vtkWriter_
->
addVertexData
(
dofvector_
,
Dune
::
VTK
::
FieldInfo
(
name_
,
VTKFieldType
<
Range
>
,
VTKFieldSize
<
Range
>
));
}
/// Implements \ref FileWriterInterface::writeFiles
virtual
void
writeFiles
(
AdaptInfo
&
adaptInfo
,
bool
force
)
override
{
using
Tree
=
typename
GlobalBasis
::
LocalView
::
Tree
;
using
Node
=
typename
Dune
::
TypeTree
::
ChildForTreePath
<
Tree
,
TreePath
>
;
writeVertexData
(
typename
Node
::
NodeTag
{},
index_
<
Node
::
CHILDREN
>
,
[
&
,
this
]()
{
if
(
mode_
==
0
)
vtkSeqWriter_
->
write
(
adaptInfo
.
getTime
(),
Dune
::
VTK
::
ascii
);
else
if
(
mode_
==
1
)
vtkSeqWriter_
->
write
(
adaptInfo
.
getTime
(),
Dune
::
VTK
::
appendedraw
);
});
vtkSeqWriter_
->
write
(
adaptInfo
.
getTime
(),
mode_
);
}
protected:
template
<
class
W
>
void
writeVertexData
(
Dune
::
TypeTree
::
LeafNodeTag
,
index_t
<
0
>
,
W
write
)
GridView
const
&
gridView
()
const
{
using
Dune
::
Functions
::
BasisBuilder
::
makeBasis
;
using
Dune
::
Functions
::
BasisBuilder
::
lagrange
;
auto
fct
=
makeDiscreteFunctionPtr
(
basis_
,
treePath_
,
vector_
);
auto
p1basis
=
makeBasis
(
basis_
->
gridView
(),
lagrange
<
1
>
());
auto
data
=
makeDOFVector
(
p1basis
,
name_
);
interpolate
(
p1basis
,
data
,
*
fct
);
auto
dataFct
=
Dune
::
Functions
::
makeDiscreteGlobalBasisFunction
<
double
>
(
p1basis
,
data
);
vtkWriter_
->
addVertexData
(
dataFct
,
Dune
::
VTK
::
FieldInfo
(
name_
,
Dune
::
VTK
::
FieldInfo
::
Type
::
scalar
,
1
));
write
();
vtkWriter_
->
clear
();
}
template
<
std
::
size_t
C
,
class
W
>
void
writeVertexData
(
Dune
::
TypeTree
::
PowerNodeTag
,
index_t
<
C
>
,
W
write
)
{
using
Dune
::
Functions
::
BasisBuilder
::
makeBasis
;
using
Dune
::
Functions
::
BasisBuilder
::
lagrange
;
using
Dune
::
Functions
::
BasisBuilder
::
power
;
using
Dune
::
Functions
::
BasisBuilder
::
flatLexicographic
;
assert
(
C
==
dow
);
auto
fct
=
makeDiscreteFunctionPtr
(
basis_
,
treePath_
,
vector_
);
auto
p1basis
=
makeBasis
(
basis_
->
gridView
(),
power
<
C
>
(
lagrange
<
1
>
(),
flatLexicographic
()));
auto
data
=
makeDOFVector
(
p1basis
,
name_
);
interpolate
(
p1basis
,
data
,
*
fct
);
using
Range
=
Dune
::
FieldVector
<
double
,
C
>
;
auto
dataFct
=
Dune
::
Functions
::
makeDiscreteGlobalBasisFunction
<
Range
>
(
p1basis
,
data
);
vtkWriter_
->
addVertexData
(
dataFct
,
Dune
::
VTK
::
FieldInfo
(
name_
,
Dune
::
VTK
::
FieldInfo
::
Type
::
vector
,
C
));
write
();
vtkWriter_
->
clear
();
return
dofvector_
.
basis
().
gridView
();
}
template
<
class
NodeTag
,
std
::
size_t
C
,
class
W
>
void
writeVertexData
(
NodeTag
,
index_t
<
C
>
,
W
)
{}
private:
std
::
shared_ptr
<
GlobalBasis
>
basis_
;
TreePath
treePath_
;
std
::
shared_ptr
<
Vector
>
vector_
;
Vector
dofvector_
;
std
::
shared_ptr
<
Dune
::
VTKWriter
<
GridView
>>
vtkWriter_
;
std
::
shared_ptr
<
Dune
::
VTKSequenceWriter
<
GridView
>>
vtkSeqWriter_
;
// std::vector<double> data_vector;
// represents VTK::OutputType: 0...ascii, 1...appendedraw
int
mode_
;
// represents VTK::OutputType: ascii, appendedraw
Dune
::
VTK
::
OutputType
mode_
;
};
template
<
class
Traits
,
class
GlobalBasis
,
class
TreePath
,
class
Vector
>
template
<
class
Traits
,
class
GlobalBasis
,
class
TreePath
>
std
::
shared_ptr
<
FileWriter
<
Traits
,
TreePath
>>
makeFileWriterPtr
(
std
::
string
baseName
,
std
::
shared_ptr
<
GlobalBasis
>
const
&
basis
,
TreePath
const
&
tp
,
std
::
shared_ptr
<
Vector
>
const
&
vector
)
DOFVectorConstView
<
GlobalBasis
,
TreePath
>
const
&
dofvector
)
{
return
std
::
make_shared
<
FileWriter
<
Traits
,
TreePath
>>
(
baseName
,
basis
,
tp
,
vector
);
return
std
::
make_shared
<
FileWriter
<
Traits
,
TreePath
>>
(
baseName
,
dof
vector
);
}
}
// end namespace AMDiS
dune/amdis/GridFunctionOperator.hpp
View file @
fa24f5f8
...
...
@@ -173,29 +173,42 @@ namespace AMDiS
:
GridFunctionOperatorBase
<
GridFct
,
QuadratureCreator
>
(
gridFct
,
0
,
quadCreator
)
{}
/// Assemble a local element vector on the element that is bound.
/// \brief Assemble a local element vector on the element that is bound.
/**
* \param contextGeometry Container for context related data
* \param quad A quadrature formula
* \param elementVector The output element vector
* \param node The node of the test-basis
**/
template
<
class
Context
,
class
QuadratureRule
,
class
ElementVector
,
class
Node
>
void
calculateElementVector
(
Context
const
&
,
//< container for context related data
QuadratureRule
const
&
,
//< a quadrature formula
ElementVector
&
,
//< the output element vector
Node
const
&
//< the node of the test-basis
)
void
calculateElementVector
(
Context
const
&
contextGeometry
,
QuadratureRule
const
&
quad
,
ElementVector
&
elementVector
,
Node
const
&
node
)
{
error_exit
(
"Needs to be implemented!"
);
}
/// Assemble a local element matrix on the element that is bound.
/// \brief Assemble a local element matrix on the element that is bound.
/**
* \param contextGeometry Container for context related data
* \param quad A quadrature formula
* \param elementMatrix The output element matrix
* \param rowNode The node of the test-basis
* \param colNode The node of the trial-basis
* \param flag1 finiteelement is the same for test and trial
* \param flag2 nodes are the same in the basis-tree
**/
template
<
class
Context
,
class
QuadratureRule
,
class
ElementMatrix
,
class
RowNode
,
class
ColNode
,
bool
sameFE
,
bool
sameNode
>
void
calculateElementMatrix
(
Context
const
&
,
//< container for context related data
QuadratureRule
const
&
,
//< a quadrature formula
ElementMatrix
&
,
//< the output element matrix
RowNode
const
&
,
//< the node of the test-basis
ColNode
const
&
,
//< the node of the trial-basis
std
::
integral_constant
<
bool
,
sameFE
>
,
//< finiteelement is the same for test and trial
std
::
integral_constant
<
bool
,
sameNode
>
//< nodes are the same in the basis-tree
)
void
calculateElementMatrix
(
Context
const
&
contextGeometry
,
QuadratureRule
const
&
quad
,
ElementMatrix
&
elementMatrix
,
RowNode
const
&
rowNode
,
ColNode
const
&
colNode
,
std
::
integral_constant
<
bool
,
sameFE
>
flag1
,
std
::
integral_constant
<
bool
,
sameNode
>
flag2
)
{
error_exit
(
"Needs to be implemented!"
);
}
...
...
@@ -261,28 +274,28 @@ namespace AMDiS
/// Store tag and expression to create a \ref GridFunctionOperator
template
<
class
Tag
,
class
Expr
>
auto
makeOperator
(
Tag
t
,
Expr
const
&
expr
)
auto
makeOperator
(
Tag
t
ag
,
Expr
const
&
expr
)
{
using
RawExpr
=
Underlying_t
<
Expr
>
;
static_assert
(
Concepts
::
HasGridFunctionOrder
<
RawExpr
>
,
"Polynomial degree of expression can not be deduced. You need to provide an explicit value for polynomial order or a quadrature rule in `makeOperator()`."
);
return
ExpressionPreOperator
<
Tag
,
Expr
,
tag
::
deduce
>
{
t
,
expr
};
return
ExpressionPreOperator
<
Tag
,
Expr
,
tag
::
deduce
>
{
t
ag
,
expr
};
}
/// Store tag and expression and polynomial order of expression to create a \ref GridFunctionOperator
template
<
class
Tag
,
class
Expr
>
auto
makeOperator
(
Tag
t
,
Expr
const
&
expr
,
int
order
,
auto
makeOperator
(
Tag
t
ag
,
Expr
const
&
expr
,
int
order
,
Dune
::
QuadratureType
::
Enum
qt
=
Dune
::
QuadratureType
::
GaussLegendre
)
{
return
ExpressionPreOperator
<
Tag
,
Expr
,
tag
::
order
>
{
t
,
expr
,
order
,
qt
};
return
ExpressionPreOperator
<
Tag
,
Expr
,
tag
::
order
>
{
t
ag
,
expr
,
order
,
qt
};
}
/// Store tag and expression and a quadrature rule to create a \ref GridFunctionOperator
template
<
class
Tag
,
class
Expr
,
class
ctype
,
int
dim
>
auto
makeOperator
(
Tag
t
,
Expr
const
&
expr
,
Dune
::
QuadratureRule
<
ctype
,
dim
>
const
&
rule
)
auto
makeOperator
(
Tag
t
ag
,
Expr
const
&
expr
,
Dune
::
QuadratureRule
<
ctype
,
dim
>
const
&
rule
)
{
return
ExpressionPreOperator
<
Tag
,
Expr
,
tag
::
rule
,
Dune
::
QuadratureRule
<
ctype
,
dim
>>
{
t
,
expr
,
rule
};
return
ExpressionPreOperator
<
Tag
,
Expr
,
tag
::
rule
,
Dune
::
QuadratureRule
<
ctype
,
dim
>>
{
t
ag
,
expr
,
rule
};
}
/** @} **/
...
...
dune/amdis/Operators.hpp
View file @
fa24f5f8
...
...
@@ -20,17 +20,17 @@
* \brief Defines operators to be assembled in the matrix/vector
*
* An `Operator` is a class providing methods necessary for assembling:
* - `
void
bind(Element, Geometry)` and `
void
unbind()` for binding an unbinding the
* - `bind(Element, Geometry)` and `unbind()` for binding an unbinding the
* element to (from) an GridView entity of codim 0. Additionally the Geometry
* object of the element is provided.
* - `
Dune::QuadratureRule<ctype,dim>
getQuadratureRule(Nodes...)` factory for the
* - `getQuadratureRule(Nodes...)` factory for the
* quadrature rules used in assembling the operator on the element. `Nodes...`
* is either `{RowNode, ColNode}` for Matrix-Operators or `{Node}` for a
* Vector-Operator.
* - `
void
calculateElementVector(ContextGeometry, QuadratureRule, ElementVector, Node)`
* - `calculateElementVector(ContextGeometry, QuadratureRule, ElementVector, Node)`
* where the `ContextGeometry` provides a reference to the ElementGeometry and
* geometry of the LocalContext (that can be different), *or*
* - `
void
calculateElementMatrix(ContextGeometry, QuadratureRule, ElementMatrix, RowNode, ColNode, Flags...)`
* - `calculateElementMatrix(ContextGeometry, QuadratureRule, ElementMatrix, RowNode, ColNode, Flags...)`
* Same as for `calculateElementVector` but additionally two optimization flags
* are provided as `bool_t<...>` type:
* + `sameFE`: the FiniteElementSpace of `RowNode` and `ColNode` are the same.
...
...
dune/amdis/ProblemStat.inc.hpp
View file @
fa24f5f8
...
...
@@ -122,8 +122,8 @@ void ProblemStat<Traits>::createFileWriter()
if
(
!
Parameters
::
get
<
std
::
string
>
(
componentName
+
"->filename"
))
return
;
auto
writer
=
makeFileWriterPtr
<
Traits
>
(
componentName
,
globalBasis
,
treePath
,
solution
);
filewriter
.
push_back
(
writer
);
auto
writer
=
makeFileWriterPtr
<
Traits
>
(
componentName
,
this
->
getSolution
(
treePath
)
);
filewriter
.
push_back
(
std
::
move
(
writer
)
)
;
});
}
...
...
dune/amdis/common/FieldMatVec.hpp
View file @
fa24f5f8
...
...
@@ -13,45 +13,48 @@
namespace
AMDiS
{
// some arithmetic operations with Dune::FieldVector
using
Dune
::
FieldVector
;
using
Dune
::
FieldMatrix
;
// some arithmetic operations with FieldVector
template
<
class
T
,
int
N
,
class
S
,
REQUIRES
(
Concepts
::
Arithmetic
<
S
>)
>
Dune
::
FieldVector
<
T
,
N
>
operator
*
(
Dune
::
FieldVector
<
T
,
N
>
v
,
S
factor
)
FieldVector
<
T
,
N
>
operator
*
(
FieldVector
<
T
,
N
>
v
,
S
factor
)
{
return
v
*=
factor
;
}
template
<
class
S
,
class
T
,
int
N
,
REQUIRES
(
Concepts
::
Arithmetic
<
S
>)
>
Dune
::
FieldVector
<
T
,
N
>
operator
*
(
S
factor
,
Dune
::
FieldVector
<
T
,
N
>
v
)
FieldVector
<
T
,
N
>
operator
*
(
S
factor
,
FieldVector
<
T
,
N
>
v
)
{
return
v
*=
factor
;
}
template
<
class
T
,
int
N
,
class
S
,
REQUIRES
(
Concepts
::
Arithmetic
<
S
>)
>
Dune
::
FieldVector
<
T
,
N
>
operator
/
(
Dune
::
FieldVector
<
T
,
N
>
v
,
S
factor
)
FieldVector
<
T
,
N
>
operator
/
(
FieldVector
<
T
,
N
>
v
,
S
factor
)
{
return
v
/=
factor
;
}
// some arithmetic operations with
Dune::
FieldMatrix
// some arithmetic operations with FieldMatrix
template
<
class
T
,
int
N
,
int
M
>
Dune
::
FieldMatrix
<
T
,
N
,
M
>
operator
+
(
Dune
::
FieldMatrix
<
T
,
N
,
M
>
A
,
Dune
::
FieldMatrix
<
T
,
N
,
M
>
const
&
B
)
FieldMatrix
<
T
,
N
,
M
>
operator
+
(
FieldMatrix
<
T
,
N
,
M
>
A
,
FieldMatrix
<
T
,
N
,
M
>
const
&
B
)
{
return
A
+=
B
;
}
template
<
class
T
,
int
N
,
int
M
>
Dune
::
FieldMatrix
<
T
,
N
,
M
>
operator
-
(
Dune
::
FieldMatrix
<
T
,
N
,
M
>
A
,
Dune
::
FieldMatrix
<
T
,
N
,
M
>
const
&
B
)
FieldMatrix
<
T
,
N
,
M
>
operator
-
(
FieldMatrix
<
T
,
N
,
M
>
A
,
FieldMatrix
<
T
,
N
,
M
>
const
&
B
)
{
return
A
-=
B
;
}
template
<
class
T
,
int
N
,
int
M
>
Dune
::
FieldVector
<
T
,
N
>
operator
*
(
Dune
::
FieldMatrix
<
T
,
N
,
M
>
const
&
mat
,
Dune
::
FieldVector
<
T
,
M
>
const
&
vec
)
FieldVector
<
T
,
N
>
operator
*
(
FieldMatrix
<
T
,
N
,
M
>
const
&
mat
,
FieldVector
<
T
,
M
>
const
&
vec
)
{
return
Dune
::
FMatrixHelp
::
mult
(
mat
,
vec
);
}
...
...
@@ -60,14 +63,14 @@ namespace AMDiS
/// Cross-product a 2d-vector = orthogonal vector
template
<
class
T
>
Dune
::
FieldVector
<
T
,
2
>
cross
(
Dune
::
FieldVector
<
T
,
2
>
const
&
a
)
FieldVector
<
T
,
2
>
cross
(
FieldVector
<
T
,
2
>
const
&
a
)
{
return
{{
a
[
1
],
-
a
[
0
]
}};
}
/// Cross-product of two vectors (in 3d only)
template
<
class
T
>
Dune
::
FieldVector
<
T
,
3
>
cross
(
Dune
::
FieldVector
<
T
,
3
>
const
&
a
,
Dune
::
FieldVector
<
T
,
3
>
const
&
b
)
FieldVector
<
T
,
3
>
cross
(
FieldVector
<
T
,
3
>
const
&
a
,
FieldVector
<
T
,
3
>
const
&
b
)
{
return
{{
a
[
1
]
*
b
[
2
]
-
a
[
2
]
*
b
[
1
],
a
[
2
]
*
b
[
0
]
-
a
[
0
]
*
b
[
2
],
...
...
@@ -76,7 +79,7 @@ namespace AMDiS
/// Dot product (vec1^T * vec2)
template
<
class
T
,
class
S
,
int
N
>
auto
dot
(
Dune
::
FieldVector
<
T
,
N
>
const
&
vec1
,
Dune
::
FieldVector
<
S
,
N
>
const
&
vec2
)
auto
dot
(
FieldVector
<
T
,
N
>
const
&
vec1
,
FieldVector
<
S
,
N
>
const
&
vec2
)
{
return
vec1
.
dot
(
vec2
);
}
...
...
@@ -86,7 +89,7 @@ namespace AMDiS
namespace
Impl
{
template
<
class
T
,
int
N
,
class
Operation
>
T
accumulate
(
Dune
::
FieldVector
<
T
,
N
>
const
&
x
,
Operation
op
)
T
accumulate
(
FieldVector
<
T
,
N
>
const
&
x
,
Operation
op
)
{
T
result
=
0
;
for
(
int
i
=
0
;
i
<
N
;
++
i
)
...
...
@@ -98,7 +101,7 @@ namespace AMDiS
/// Sum of vector entires.
template
<
class
T
,
int
N
>
T
sum
(
Dune
::
FieldVector
<
T
,
N
>
const
&
x
)
T
sum
(
FieldVector
<
T
,
N
>
const
&
x
)
{
return
Impl
::
accumulate
(
x
,
Operation
::
Plus
{});
}
...
...
@@ -106,7 +109,7 @@ namespace AMDiS
/// Dot-product with the vector itself
template
<
class
T
,
int
N
>
auto
unary_dot
(
Dune
::
FieldVector
<
T
,
N
>
const
&
x
)
auto
unary_dot
(
FieldVector
<
T
,
N
>
const
&
x
)
{
auto
op
=
[](
auto
const
&
a
,
auto
const
&
b
)
{
return
a
+
Math
::
sqr
(
std
::
abs
(
b
));
};
return
Impl
::
accumulate
(
x
,
op
);
...
...
@@ -114,28 +117,28 @@ namespace AMDiS
/// Maximum over all vector entries
template
<
class
T
,
int
N
>
auto
max
(
Dune
::
FieldVector
<
T
,
N
>
const
&
x
)
auto
max
(
FieldVector
<
T
,
N
>
const
&
x
)
{
return
Impl
::
accumulate
(
x
,
Operation
::
Max
{});
}
/// Minimum over all vector entries
template
<
class
T
,
int
N
>
auto
min
(
Dune
::
FieldVector
<
T
,
N
>
const
&
x
)
auto
min
(
FieldVector
<
T
,
N
>
const
&
x
)
{
return
Impl
::
accumulate
(
x
,
Operation
::
Min
{});
}
/// Maximum of the absolute values of vector entries
template
<
class
T
,
int
N
>
auto
abs_max
(
Dune
::
FieldVector
<
T
,
N
>
const
&
x
)
auto
abs_max
(
FieldVector
<
T
,
N
>
const
&
x
)
{
return
Impl
::
accumulate
(
x
,
Operation
::
AbsMax
{});
}
/// Minimum of the absolute values of vector entries
template
<
class
T
,
int
N
>
auto
abs_min
(
Dune
::
FieldVector
<
T
,
N
>
const
&
x
)
auto
abs_min
(
FieldVector
<
T
,
N
>
const
&
x
)
{
return
Impl
::
accumulate
(
x
,
Operation
::
AbsMin
{});
}
...
...
@@ -146,7 +149,7 @@ namespace AMDiS
* \brief The 1-norm of a vector = sum_i |x_i|
**/
template
<
class
T
,
int
N
>
auto
one_norm
(
Dune
::
FieldVector
<
T
,
N
>
const
&
x
)
auto
one_norm
(
FieldVector
<
T
,
N
>
const
&
x
)
{
auto
op
=
[](
auto
const
&
a
,
auto
const
&
b
)
{
return
a
+
std
::
abs
(
b
);
};
return
Impl
::
accumulate
(
x
,
op
);
...
...
@@ -156,7 +159,7 @@ namespace AMDiS
* \brief The euklidean 2-norm of a vector = sqrt( sum_i |x_i|^2 )
**/
template
<
class
T
,
int
N
>
auto
two_norm
(
Dune
::
FieldVector
<
T
,
N
>
const
&
x
)
auto
two_norm
(
FieldVector
<
T
,
N
>
const
&
x
)
{
return
std
::
sqrt
(
unary_dot
(
x
));
}
...
...
@@ -165,7 +168,7 @@ namespace AMDiS
* \brief The p-norm of a vector = ( sum_i |x_i|^p )^(1/p)
**/
template
<
int
p
,
class
T
,
int
N
>
auto
p_norm
(
Dune
::
FieldVector
<
T
,
N
>
const
&
x
)
auto
p_norm
(
FieldVector
<
T
,
N
>
const
&
x
)
{
auto
op
=
[](
auto
const
&
a
,
auto
const
&
b
)
{
return
a
+
Math
::
pow
<
p
>
(
std
::
abs
(
b
));
};
return
std
::
pow
(
Impl
::
accumulate
(
x
,
op
),
1.0
/
p
);
...
...
@@ -175,7 +178,7 @@ namespace AMDiS
* \brief The infty-norm of a vector = max_i |x_i| = alias for \ref abs_max
**/
template
<
class
T
,
int
N
>
auto
infty_norm
(
Dune
::
FieldVector
<
T
,
N
>
const
&
x
)
auto
infty_norm
(
FieldVector
<
T
,
N
>
const
&
x
)
{
return
abs_max
(
x
);
}
...
...
@@ -184,7 +187,7 @@ namespace AMDiS
/// The euklidean distance between two vectors = |lhs-rhs|_2
template
<
class
T
,
int
N
>
T
distance
(
Dune
::
FieldVector
<
T
,
N
>
const
&
lhs
,
Dune
::
FieldVector
<
T
,
N
>
const
&
rhs
)
T
distance
(
FieldVector
<
T
,
N
>
const
&
lhs
,
FieldVector
<
T
,
N
>
const
&
rhs
)
{
T
result
=
0
;
for
(
int
i
=
0
;
i
<
N
;
++
i
)
...
...
@@ -196,9 +199,9 @@ namespace AMDiS
/// Outer product (vec1 * vec2^T)
template
<
class
T
,
class
S
,
int
N
,
int
M
,
int
K
>
auto
outer
(
Dune
::
FieldMatrix
<
T
,
N
,
K
>
const
&
vec1
,
Dune
::
FieldMatrix
<
S
,
M
,
K
>
const
&
vec2
)
auto
outer
(
FieldMatrix
<
T
,
N
,
K
>
const
&
vec1
,
FieldMatrix
<
S
,
M
,
K
>
const
&
vec2
)
{
using
result_type
=
Dune
::
FieldMatrix
<
decltype
(
std
::
declval
<
T
>
()
*
std
::
declval
<
S
>
()
),
N
,
M
>
;
using
result_type
=
FieldMatrix
<
decltype
(
std
::
declval
<
T
>
()
*
std
::
declval
<
S
>
()
),
N
,
M
>
;
result_type
mat
;
for
(
int
i
=
0
;
i
<
N
;
++
i
)
for
(
int
j
=
0
;
j
<
M
;
++
j
)
...
...
@@ -209,28 +212,28 @@ namespace AMDiS
// ----------------------------------------------------------------------------
template
<
class
T
>
T
det
(
Dune
::
FieldMatrix
<
T
,
0
,
0
>
const
&
/*mat*/
)
T
det
(
FieldMatrix
<
T
,
0
,
0
>
const
&
/*mat*/
)
{
return
0
;
}
/// Determinant of a 1x1 matrix
template
<
class
T
>
T
det
(
Dune
::
FieldMatrix
<
T
,
1
,
1
>
const
&
mat
)
T
det
(
FieldMatrix
<
T
,
1
,
1
>
const
&
mat
)
{
return
mat
[
0
][
0
];
}
/// Determinant of a 2x2 matrix
template
<
class
T
>
T
det
(
Dune
::
FieldMatrix
<
T
,
2
,
2
>
const
&
mat
)
T
det
(
FieldMatrix
<
T
,
2
,
2
>
const
&
mat
)
{
return
mat
[
0
][
0
]
*
mat
[
1
][
1
]
-
mat
[
0
][
1
]
*
mat
[
1
][
0
];
}
/// Determinant of a 3x3 matrix
template
<
class
T
>
T
det
(
Dune
::
FieldMatrix
<
T
,
3
,
3
>
const
&
mat
)
T
det
(
FieldMatrix
<
T
,
3
,
3
>
const
&
mat
)
{
return
mat
[
0
][
0
]
*
mat
[
1
][
1
]
*
mat
[
2
][
2
]
+
mat
[
0
][
1
]
*
mat
[
1
][
2
]
*
mat
[
2
][
0
]
+
mat
[
0
][
2
]
*
mat
[
1
][
0
]
*
mat
[
2
][
1
]
-
(
mat
[
0
][
2
]
*
mat
[
1
][
1
]
*
mat
[
2
][
0
]
+
mat
[
0
][
1
]
*
mat
[
1
][
0
]
*
mat
[
2
][
2
]
+
mat
[
0
][
0
]
*
mat
[
1
][
2
]
*
mat
[
2
][
1
]);
...
...
@@ -238,14 +241,14 @@ namespace AMDiS
/// Determinant of a NxN matrix
template
<
class
T
,
int
N
>
T
det
(
Dune
::
FieldMatrix
<
T
,
N
,
N
>
const
&
mat
)
T
det
(
FieldMatrix
<
T
,
N
,
N
>
const
&
mat
)
{
return
mat
.
determinant
();