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
3d0fdf35
Commit
3d0fdf35
authored
Dec 04, 2018
by
Praetorius, Simon
Browse files
added some missing operators and a test for all operators
parent
746bb31d
Changes
8
Hide whitespace changes
Inline
Side-by-side
src/amdis/Operators.hpp
View file @
3d0fdf35
...
...
@@ -46,6 +46,10 @@
#include
<amdis/assembler/ZeroOrderTestvecTrialvec.hpp>
// <Psi, A * Phi>, <Psi, c * Phi>
// first-order operators
#include
<amdis/assembler/FirstOrderPartialTest.hpp>
// <d_i(psi), c>
#include
<amdis/assembler/FirstOrderGradTest.hpp>
// <grad(psi), b>
#include
<amdis/assembler/FirstOrderDivTestvec.hpp>
// <div(Psi), c>
#include
<amdis/assembler/FirstOrderDivTestvecTrial.hpp>
// <div(Psi), c * phi>
#include
<amdis/assembler/FirstOrderGradTestTrial.hpp>
// <grad(psi), b * phi>
#include
<amdis/assembler/FirstOrderGradTestTrialvec.hpp>
// <grad(psi), c * Phi>
...
...
src/amdis/assembler/ConvectionDiffusionOperator.hpp
View file @
3d0fdf35
...
...
@@ -108,11 +108,11 @@ namespace AMDiS
const
auto
c
=
localFctC
(
local
);
IF_CONSTEXPR
(
conserving
)
{
WorldVector
gradAi
,
gradBi
;
WorldVector
gradAi
;
for
(
std
::
size_t
i
=
0
;
i
<
numLocalFe
;
++
i
)
{
const
int
local_i
=
rowNode
.
localIndex
(
i
);
gradAi
=
A
*
gradients
[
i
];
gradBi
=
b
*
gradients
[
i
];
auto
gradBi
=
b
*
gradients
[
i
];
for
(
std
::
size_t
j
=
0
;
j
<
numLocalFe
;
++
j
)
{
const
int
local_j
=
colNode
.
localIndex
(
j
);
elementMatrix
[
local_i
][
local_j
]
+=
(
dot
(
gradAi
,
gradients
[
j
])
+
(
c
*
shapeValues
[
i
]
-
gradBi
)
*
shapeValues
[
j
])
*
factor
;
...
...
src/amdis/assembler/FirstOrderDivTestvec.hpp
0 → 100644
View file @
3d0fdf35
#pragma once
#include
<type_traits>
#include
<amdis/GridFunctionOperator.hpp>
#include
<amdis/LocalBasisCache.hpp>
#include
<amdis/Output.hpp>
#include
<amdis/common/Size.hpp>
namespace
AMDiS
{
/**
* \addtogroup operators
* @{
**/
namespace
tag
{
struct
divtestvec
{};
}
/// first-order operator \f$ \langle\nabla\cdot\Psi, c\rangle \f$
template
<
class
LocalContext
,
class
GridFct
>
class
GridFunctionOperator
<
tag
::
divtestvec
,
LocalContext
,
GridFct
>
:
public
GridFunctionOperatorBase
<
GridFunctionOperator
<
tag
::
divtestvec
,
LocalContext
,
GridFct
>
,
LocalContext
,
GridFct
>
{
using
Self
=
GridFunctionOperator
;
using
Super
=
GridFunctionOperatorBase
<
Self
,
LocalContext
,
GridFct
>
;
static_assert
(
Size_v
<
typename
GridFct
::
Range
>
==
1
,
"Expression must be of scalar type."
);
public:
GridFunctionOperator
(
tag
::
divtestvec
,
GridFct
const
&
expr
)
:
Super
(
expr
,
1
)
{}
template
<
class
Context
,
class
Node
,
class
ElementVector
>
void
getElementVector
(
Context
const
&
context
,
Node
const
&
node
,
ElementVector
&
elementVector
)
{
static_assert
(
Node
::
isPower
,
"Node must be a Power-Node."
);
static
const
std
::
size_t
CHILDREN
=
Node
::
CHILDREN
;
static_assert
(
CHILDREN
==
Context
::
dow
,
""
);
auto
const
&
quad
=
this
->
getQuadratureRule
(
context
.
type
(),
node
);
auto
const
&
localFE
=
node
.
child
(
0
).
finiteElement
();
std
::
size_t
feSize
=
localFE
.
size
();
using
LocalBasisType
=
typename
std
::
decay_t
<
decltype
(
localFE
)
>::
Traits
::
LocalBasisType
;
using
RangeFieldType
=
typename
LocalBasisType
::
Traits
::
RangeFieldType
;
LocalBasisCache
<
LocalBasisType
>
cache
(
localFE
.
localBasis
());
auto
const
&
shapeGradientsCache
=
cache
.
evaluateJacobianAtQp
(
context
,
quad
);
for
(
std
::
size_t
iq
=
0
;
iq
<
quad
.
size
();
++
iq
)
{
// Position of the current quadrature point in the reference element
decltype
(
auto
)
local
=
context
.
local
(
quad
[
iq
].
position
());
// The transposed inverse Jacobian of the map from the reference element to the element
const
auto
jacobian
=
context
.
geometry
().
jacobianInverseTransposed
(
local
);
// The multiplicative factor in the integral transformation formula
const
auto
factor
=
Super
::
coefficient
(
local
)
*
context
.
integrationElement
(
quad
[
iq
].
position
())
*
quad
[
iq
].
weight
();
// The gradients of the shape functions on the reference element
auto
const
&
shapeGradients
=
shapeGradientsCache
[
iq
];
// Compute the shape function gradients on the real element
using
WorldVector
=
FieldVector
<
RangeFieldType
,
Context
::
dow
>
;
thread_local
std
::
vector
<
WorldVector
>
gradients
;
gradients
.
resize
(
shapeGradients
.
size
());
for
(
std
::
size_t
i
=
0
;
i
<
gradients
.
size
();
++
i
)
jacobian
.
mv
(
shapeGradients
[
i
][
0
],
gradients
[
i
]);
for
(
std
::
size_t
j
=
0
;
j
<
feSize
;
++
j
)
{
for
(
std
::
size_t
k
=
0
;
k
<
CHILDREN
;
++
k
)
{
const
auto
local_kj
=
node
.
child
(
k
).
localIndex
(
j
);
elementVector
[
local_kj
]
+=
factor
*
gradients
[
j
][
k
];
}
}
}
}
};
/** @} **/
}
// end namespace AMDiS
src/amdis/assembler/FirstOrderGradTest.hpp
0 → 100644
View file @
3d0fdf35
#pragma once
#include
<type_traits>
#include
<amdis/GridFunctionOperator.hpp>
#include
<amdis/LocalBasisCache.hpp>
#include
<amdis/Output.hpp>
#include
<amdis/common/Size.hpp>
namespace
AMDiS
{
/**
* \addtogroup operators
* @{
**/
namespace
tag
{
struct
gradtest
{};
}
/// first-order operator \f$ \langle\Psi, c\,\nabla\phi\rangle \f$
template
<
class
LocalContext
,
class
GridFct
>
class
GridFunctionOperator
<
tag
::
gradtest
,
LocalContext
,
GridFct
>
:
public
GridFunctionOperatorBase
<
GridFunctionOperator
<
tag
::
gradtest
,
LocalContext
,
GridFct
>
,
LocalContext
,
GridFct
>
{
static
const
int
dow
=
ContextGeometry
<
LocalContext
>::
dow
;
using
Self
=
GridFunctionOperator
;
using
Super
=
GridFunctionOperatorBase
<
Self
,
LocalContext
,
GridFct
>
;
static_assert
(
Size_v
<
typename
GridFct
::
Range
>
==
dow
,
"Expression must be of vector type."
);
public:
GridFunctionOperator
(
tag
::
gradtest
,
GridFct
const
&
expr
)
:
Super
(
expr
,
1
)
{}
template
<
class
Context
,
class
Node
,
class
ElementVector
>
void
getElementVector
(
Context
const
&
context
,
Node
const
&
node
,
ElementVector
&
elementVector
)
{
static_assert
(
Node
::
isLeaf
,
"Node must be Leaf-Node."
);
auto
const
&
quad
=
this
->
getQuadratureRule
(
context
.
type
(),
node
);
auto
const
&
localFE
=
node
.
finiteElement
();
std
::
size_t
feSize
=
localFE
.
size
();
using
LocalBasisType
=
typename
std
::
decay_t
<
decltype
(
localFE
)
>::
Traits
::
LocalBasisType
;
using
RangeFieldType
=
typename
LocalBasisType
::
Traits
::
RangeFieldType
;
LocalBasisCache
<
LocalBasisType
>
cache
(
localFE
.
localBasis
());
auto
const
&
shapeGradientsCache
=
cache
.
evaluateJacobianAtQp
(
context
,
quad
);
for
(
std
::
size_t
iq
=
0
;
iq
<
quad
.
size
();
++
iq
)
{
// Position of the current quadrature point in the reference element
decltype
(
auto
)
local
=
context
.
local
(
quad
[
iq
].
position
());
// The transposed inverse Jacobian of the map from the reference element to the element
const
auto
jacobian
=
context
.
geometry
().
jacobianInverseTransposed
(
local
);
// The multiplicative factor in the integral transformation formula
const
auto
factor
=
Super
::
coefficient
(
local
);
const
auto
dx
=
context
.
integrationElement
(
quad
[
iq
].
position
())
*
quad
[
iq
].
weight
();
// The gradients of the shape functions on the reference element
auto
const
&
shapeGradients
=
shapeGradientsCache
[
iq
];
// Compute the shape function gradients on the real element
using
WorldVector
=
FieldVector
<
RangeFieldType
,
Context
::
dow
>
;
thread_local
std
::
vector
<
WorldVector
>
gradients
;
gradients
.
resize
(
shapeGradients
.
size
());
for
(
std
::
size_t
i
=
0
;
i
<
gradients
.
size
();
++
i
)
jacobian
.
mv
(
shapeGradients
[
i
][
0
],
gradients
[
i
]);
for
(
std
::
size_t
i
=
0
;
i
<
feSize
;
++
i
)
{
const
auto
local_i
=
node
.
localIndex
(
i
);
elementVector
[
local_i
]
+=
dx
*
(
factor
*
gradients
[
i
]);
}
}
}
};
/** @} **/
}
// end namespace AMDiS
src/amdis/assembler/FirstOrderPartialTest.hpp
0 → 100644
View file @
3d0fdf35
#pragma once
#include
<type_traits>
#include
<amdis/assembler/FirstOrderTestPartialTrial.hpp>
namespace
AMDiS
{
/**
* \addtogroup operators
* @{
**/
namespace
tag
{
struct
partialtest
{
std
::
size_t
comp
;
};
}
/// first-order operator \f$ \langle\partial_i\psi, c\rangle \f$
template
<
class
LocalContext
,
class
GridFct
>
class
GridFunctionOperator
<
tag
::
partialtest
,
LocalContext
,
GridFct
>
:
public
GridFunctionOperatorBase
<
GridFunctionOperator
<
tag
::
partialtest
,
LocalContext
,
GridFct
>
,
LocalContext
,
GridFct
>
{
using
Self
=
GridFunctionOperator
;
using
Super
=
GridFunctionOperatorBase
<
Self
,
LocalContext
,
GridFct
>
;
static_assert
(
Size_v
<
typename
GridFct
::
Range
>
==
1
,
"Expression must be of scalar type."
);
public:
GridFunctionOperator
(
tag
::
partialtest
tag
,
GridFct
const
&
expr
)
:
Super
(
expr
,
1
)
,
comp_
(
tag
.
comp
)
{}
template
<
class
Context
,
class
Node
,
class
ElementVector
>
void
getElementVector
(
Context
const
&
context
,
Node
const
&
node
,
ElementVector
&
elementVector
)
{
static_assert
(
Node
::
isLeaf
,
"Operator can be applied to Leaf-Nodes only."
);
auto
const
&
quad
=
this
->
getQuadratureRule
(
context
.
type
(),
node
);
auto
const
&
localFE
=
node
.
finiteElement
();
std
::
size_t
feSize
=
localFE
.
size
();
using
LocalBasisType
=
typename
std
::
decay_t
<
decltype
(
localFE
)
>::
Traits
::
LocalBasisType
;
using
RangeFieldType
=
typename
LocalBasisType
::
Traits
::
RangeFieldType
;
LocalBasisCache
<
LocalBasisType
>
cache
(
localFE
.
localBasis
());
auto
const
&
shapeGradientsCache
=
cache
.
evaluateJacobianAtQp
(
context
,
quad
);
for
(
std
::
size_t
iq
=
0
;
iq
<
quad
.
size
();
++
iq
)
{
// Position of the current quadrature point in the reference element
decltype
(
auto
)
local
=
context
.
local
(
quad
[
iq
].
position
());
// The transposed inverse Jacobian of the map from the reference element to the element
const
auto
jacobian
=
context
.
geometry
().
jacobianInverseTransposed
(
local
);
assert
(
jacobian
.
M
()
==
Context
::
dim
);
// The multiplicative factor in the integral transformation formula
const
auto
dx
=
context
.
integrationElement
(
quad
[
iq
].
position
())
*
quad
[
iq
].
weight
();
const
auto
c
=
Super
::
coefficient
(
local
);
// The gradients of the shape functions on the reference element
auto
const
&
shapeGradients
=
shapeGradientsCache
[
iq
];
// Compute the shape function gradients on the real element
thread_local
std
::
vector
<
RangeFieldType
>
partial
;
partial
.
resize
(
shapeGradients
.
size
());
for
(
std
::
size_t
i
=
0
;
i
<
partial
.
size
();
++
i
)
{
partial
[
i
]
=
jacobian
[
comp_
][
0
]
*
shapeGradients
[
i
][
0
][
0
];
for
(
std
::
size_t
j
=
1
;
j
<
jacobian
.
M
();
++
j
)
partial
[
i
]
+=
jacobian
[
comp_
][
j
]
*
shapeGradients
[
i
][
0
][
j
];
}
for
(
std
::
size_t
j
=
0
;
j
<
feSize
;
++
j
)
{
const
auto
local_j
=
node
.
localIndex
(
j
);
elementVector
[
local_j
]
+=
dx
*
(
c
*
partial
[
j
]);
}
}
}
private:
std
::
size_t
comp_
;
};
/** @} **/
}
// end namespace AMDiS
src/amdis/assembler/SecondOrderPartialTestPartialTrial.hpp
View file @
3d0fdf35
...
...
@@ -88,7 +88,7 @@ namespace AMDiS
rowPartial
.
resize
(
rowShapeGradients
.
size
());
for
(
std
::
size_t
i
=
0
;
i
<
rowPartial
.
size
();
++
i
)
{
rowPartial
[
i
]
=
jacobian
[
compTest_
][
0
]
*
rowShapeGradients
[
i
][
0
][
0
];
for
(
std
::
size_t
j
=
1
;
j
<
jacobian
.
cols
()
;
++
j
)
for
(
std
::
size_t
j
=
1
;
j
<
jacobian
.
cols
;
++
j
)
rowPartial
[
i
]
+=
jacobian
[
compTest_
][
j
]
*
rowShapeGradients
[
i
][
0
][
j
];
}
...
...
@@ -96,7 +96,7 @@ namespace AMDiS
colPartial
.
resize
(
colShapeGradients
.
size
());
for
(
std
::
size_t
i
=
0
;
i
<
colPartial
.
size
();
++
i
)
{
colPartial
[
i
]
=
jacobian
[
compTrial_
][
0
]
*
colShapeGradients
[
i
][
0
][
0
];
for
(
std
::
size_t
j
=
1
;
j
<
jacobian
.
cols
()
;
++
j
)
for
(
std
::
size_t
j
=
1
;
j
<
jacobian
.
cols
;
++
j
)
colPartial
[
i
]
+=
jacobian
[
compTrial_
][
j
]
*
colShapeGradients
[
i
][
0
][
j
];
}
...
...
test/CMakeLists.txt
View file @
3d0fdf35
...
...
@@ -40,6 +40,9 @@ dune_add_test(SOURCES MultiTypeVectorTest.cpp
dune_add_test
(
SOURCES MultiTypeMatrixTest.cpp
LINK_LIBRARIES amdis
)
dune_add_test
(
SOURCES OperatorsTest.cpp
LINK_LIBRARIES amdis
)
dune_add_test
(
SOURCES RangeTypeTest.cpp
LINK_LIBRARIES amdis
)
...
...
test/OperatorsTest.cpp
0 → 100644
View file @
3d0fdf35
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include
"config.h"
#include
<amdis/AMDiS.hpp>
#include
<amdis/ProblemStat.hpp>
#include
<amdis/Operators.hpp>
#include
<amdis/assembler/ConvectionDiffusionOperator.hpp>
#include
<amdis/assembler/StokesOperator.hpp>
using
namespace
AMDiS
;
using
Grid
=
Dune
::
YaspGrid
<
2
>
;
using
Param
=
TaylorHoodBasis
<
typename
Grid
::
LeafGridView
>
;
using
Problem
=
ProblemStat
<
Param
>
;
int
main
(
int
argc
,
char
**
argv
)
{
AMDiS
::
init
(
argc
,
argv
);
using
namespace
Dune
::
Indices
;
Problem
prob
(
"ellipt"
);
prob
.
initialize
(
INIT_ALL
);
auto
_v
=
_0
;
auto
_p
=
_1
;
FieldVector
<
double
,
1
>
scalar
=
1.0
;
FieldVector
<
double
,
2
>
vec
;
vec
=
1.0
;
FieldMatrix
<
double
,
2
,
2
>
mat
;
mat
=
1.0
;
auto
op1
=
makeOperator
(
tag
::
divtestvec
{},
scalar
);
prob
.
addVectorOperator
(
op1
,
_v
);
auto
op2
=
makeOperator
(
tag
::
divtestvec_trial
{},
scalar
);
prob
.
addMatrixOperator
(
op2
,
_v
,
_p
);
auto
op3
=
makeOperator
(
tag
::
gradtest
{},
vec
);
prob
.
addVectorOperator
(
op3
,
_p
);
auto
op4
=
makeOperator
(
tag
::
gradtest_trial
{},
vec
);
prob
.
addMatrixOperator
(
op4
,
_p
,
_p
);
auto
op5a
=
makeOperator
(
tag
::
gradtest_trialvec
{},
scalar
);
// auto op5b = makeOperator(tag::gradtest_trialvec{}, mat); // TODO: needs to be implemented
prob
.
addMatrixOperator
(
op5a
,
_p
,
_v
);
auto
op6
=
makeOperator
(
tag
::
partialtest
{
0
},
scalar
);
prob
.
addVectorOperator
(
op6
,
_p
);
auto
op7
=
makeOperator
(
tag
::
partialtest_trial
{
0
},
scalar
);
prob
.
addMatrixOperator
(
op7
,
_p
,
_p
);
auto
op8
=
makeOperator
(
tag
::
test_divtrialvec
{},
scalar
);
prob
.
addMatrixOperator
(
op8
,
_p
,
_v
);
auto
op9
=
makeOperator
(
tag
::
test_gradtrial
{},
vec
);
prob
.
addMatrixOperator
(
op9
,
_p
,
_p
);
auto
op10
=
makeOperator
(
tag
::
test_partialtrial
{
0
},
scalar
);
prob
.
addMatrixOperator
(
op10
,
_p
,
_p
);
auto
op11a
=
makeOperator
(
tag
::
testvec_gradtrial
{},
scalar
);
// auto op11b = makeOperator(tag::testvec_gradtrial{}, mat); // TODO: needs to be implemented
prob
.
addMatrixOperator
(
op11a
,
_v
,
_p
);
auto
op12
=
makeOperator
(
tag
::
divtestvec_divtrialvec
{},
scalar
);
prob
.
addMatrixOperator
(
op12
,
_v
,
_v
);
auto
op13a
=
makeOperator
(
tag
::
gradtest_gradtrial
{},
scalar
);
auto
op13b
=
makeOperator
(
tag
::
gradtest_gradtrial
{},
mat
);
prob
.
addMatrixOperator
(
op13a
,
_p
,
_p
);
prob
.
addMatrixOperator
(
op13b
,
_p
,
_p
);
auto
op14
=
makeOperator
(
tag
::
partialtest_partialtrial
{
0
,
0
},
scalar
);
prob
.
addMatrixOperator
(
op14
,
_p
,
_p
);
auto
op15
=
makeOperator
(
tag
::
test
{},
scalar
);
prob
.
addVectorOperator
(
op15
,
_p
);
auto
op16
=
makeOperator
(
tag
::
test_trial
{},
scalar
);
prob
.
addMatrixOperator
(
op16
,
_p
,
_p
);
auto
op17
=
makeOperator
(
tag
::
test_trialvec
{},
vec
);
prob
.
addMatrixOperator
(
op17
,
_p
,
_v
);
auto
op18
=
makeOperator
(
tag
::
testvec
{},
vec
);
prob
.
addVectorOperator
(
op18
,
_v
);
auto
op19
=
makeOperator
(
tag
::
testvec_trial
{},
vec
);
prob
.
addMatrixOperator
(
op19
,
_v
,
_p
);
auto
op20a
=
makeOperator
(
tag
::
testvec_trialvec
{},
scalar
);
auto
op20b
=
makeOperator
(
tag
::
testvec_trialvec
{},
mat
);
prob
.
addMatrixOperator
(
op20a
,
_v
,
_v
);
prob
.
addMatrixOperator
(
op20b
,
_v
,
_v
);
// some special operators
auto
opStokes
=
makeOperator
(
tag
::
stokes
{},
scalar
);
prob
.
addMatrixOperator
(
opStokes
);
prob
.
addMatrixOperator
(
opStokes
,
{},
{});
auto
opCDa
=
convectionDiffusion
(
scalar
,
scalar
,
scalar
,
scalar
);
prob
.
addMatrixOperator
(
opCDa
,
_p
,
_p
);
prob
.
addVectorOperator
(
opCDa
,
_p
);
auto
opCDb
=
convectionDiffusion
(
mat
,
vec
,
scalar
,
scalar
,
std
::
false_type
{});
prob
.
addMatrixOperator
(
opCDb
,
_p
,
_p
);
prob
.
addVectorOperator
(
opCDb
,
_p
);
AMDiS
::
finalize
();
return
0
;
}
Write
Preview
Supports
Markdown
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