Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
iwr
amdis
Commits
30debe8f
Commit
30debe8f
authored
Feb 24, 2012
by
Praetorius, Simon
Browse files
several changes in integrate routines of DOFVectors
parent
63c967c1
Changes
6
Hide whitespace changes
Inline
Side-by-side
AMDiS/CMakeLists.txt
View file @
30debe8f
...
...
@@ -8,8 +8,8 @@ SET(SOURCE_DIR ${AMDIS_SOURCE_DIR}/src)
#TODO: use the cmake build type
SET
(
MTL_INCLUDE_DIR
${
LIB_DIR
}
/mtl4/ CACHE PATH
"mtl4 directory"
)
set
(
CMAKE_CXX_FLAGS_RELEASE
"
${
CMAKE_CXX_FLAGS_RELEASE
}
-g -Wall -DDEBUG=0"
)
set
(
CMAKE_CXX_FLAGS_DEBUG
"
${
CMAKE_CXX_FLAGS_DEBUG
}
-Wall -DDEBUG=1 -O0"
)
set
(
CMAKE_CXX_FLAGS_RELEASE
"
${
CMAKE_CXX_FLAGS_RELEASE
}
-g -Wall
-Wno-unused-but-set-variable
-DDEBUG=0"
)
set
(
CMAKE_CXX_FLAGS_DEBUG
"
${
CMAKE_CXX_FLAGS_DEBUG
}
-Wall
-Wno-unused-but-set-variable
-DDEBUG=1 -O0"
)
if
(
NOT CurrentRevision
)
find_package
(
Subversion
)
if
(
Subversion_FOUND
)
...
...
AMDiS/src/DOFVector.cc
View file @
30debe8f
...
...
@@ -714,193 +714,6 @@ namespace AMDiS {
}
double
integrate
(
const
DOFVector
<
double
>
&
vec1
,
const
DOFVector
<
double
>
&
vec2
,
BinaryAbstractFunction
<
double
,
double
,
double
>
*
fct
)
{
if
(
vec1
.
getFeSpace
()
->
getMesh
()
==
vec2
.
getFeSpace
()
->
getMesh
())
return
intSingleMesh
(
vec1
,
vec2
,
fct
);
else
return
intDualMesh
(
vec1
,
vec2
,
fct
);
}
double
intSingleMesh
(
const
DOFVector
<
double
>
&
vec1
,
const
DOFVector
<
double
>
&
vec2
,
BinaryAbstractFunction
<
double
,
double
,
double
>
*
fct
)
{
FUNCNAME
(
"intDualmesh()"
);
TEST_EXIT
(
fct
)(
"No function defined!
\n
"
);
int
deg
=
std
::
max
(
fct
->
getDegree
(),
2
*
vec1
.
getFeSpace
()
->
getBasisFcts
()
->
getDegree
());
Quadrature
*
quad
=
Quadrature
::
provideQuadrature
(
vec1
.
getFeSpace
()
->
getMesh
()
->
getDim
(),
deg
);
FastQuadrature
*
fastQuad
=
FastQuadrature
::
provideFastQuadrature
(
vec1
.
getFeSpace
()
->
getBasisFcts
(),
*
quad
,
INIT_PHI
);
mtl
::
dense_vector
<
double
>
qp1
(
fastQuad
->
getNumPoints
());
mtl
::
dense_vector
<
double
>
qp2
(
fastQuad
->
getNumPoints
());
double
value
=
0.0
;
Flag
traverseFlag
=
Mesh
::
CALL_LEAF_EL
|
Mesh
::
FILL_COORDS
|
Mesh
::
FILL_DET
;
TraverseStack
stack
;
ElInfo
*
elInfo
=
stack
.
traverseFirst
(
vec1
.
getFeSpace
()
->
getMesh
(),
-
1
,
traverseFlag
);
while
(
elInfo
)
{
vec1
.
getVecAtQPs
(
elInfo
,
quad
,
fastQuad
,
qp1
);
vec2
.
getVecAtQPs
(
elInfo
,
quad
,
fastQuad
,
qp2
);
double
tmp
=
0.0
;
for
(
int
iq
=
0
;
iq
<
fastQuad
->
getNumPoints
();
iq
++
)
tmp
+=
fastQuad
->
getWeight
(
iq
)
*
(
*
fct
)(
qp1
[
iq
],
qp2
[
iq
]);
value
+=
tmp
*
elInfo
->
getDet
();
elInfo
=
stack
.
traverseNext
(
elInfo
);
}
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
double
localValue
=
value
;
MPI
::
COMM_WORLD
.
Allreduce
(
&
localValue
,
&
value
,
1
,
MPI_DOUBLE
,
MPI_SUM
);
#endif
return
value
;
}
double
intDualMesh
(
const
DOFVector
<
double
>
&
vec1
,
const
DOFVector
<
double
>
&
vec2
,
BinaryAbstractFunction
<
double
,
double
,
double
>
*
fct
)
{
FUNCNAME
(
"intDualmesh()"
);
TEST_EXIT
(
fct
)(
"No function defined!
\n
"
);
int
deg
=
std
::
max
(
fct
->
getDegree
(),
2
*
vec1
.
getFeSpace
()
->
getBasisFcts
()
->
getDegree
());
Quadrature
*
quad
=
Quadrature
::
provideQuadrature
(
vec1
.
getFeSpace
()
->
getMesh
()
->
getDim
(),
deg
);
FastQuadrature
*
fastQuad
=
FastQuadrature
::
provideFastQuadrature
(
vec1
.
getFeSpace
()
->
getBasisFcts
(),
*
quad
,
INIT_PHI
);
mtl
::
dense_vector
<
double
>
qp1
(
fastQuad
->
getNumPoints
());
mtl
::
dense_vector
<
double
>
qp2
(
fastQuad
->
getNumPoints
());
double
value
=
0.0
;
Flag
traverseFlag
=
Mesh
::
CALL_LEAF_EL
|
Mesh
::
FILL_COORDS
|
Mesh
::
FILL_DET
;
DualTraverse
dualTraverse
;
DualElInfo
dualElInfo
;
bool
cont
=
dualTraverse
.
traverseFirst
(
vec1
.
getFeSpace
()
->
getMesh
(),
vec2
.
getFeSpace
()
->
getMesh
(),
-
1
,
-
1
,
traverseFlag
,
traverseFlag
,
dualElInfo
);
while
(
cont
)
{
vec1
.
getVecAtQPs
(
dualElInfo
.
smallElInfo
,
dualElInfo
.
largeElInfo
,
quad
,
NULL
,
qp1
);
vec2
.
getVecAtQPs
(
dualElInfo
.
smallElInfo
,
dualElInfo
.
largeElInfo
,
quad
,
NULL
,
qp2
);
double
tmp
=
0.0
;
for
(
int
iq
=
0
;
iq
<
quad
->
getNumPoints
();
iq
++
)
tmp
+=
quad
->
getWeight
(
iq
)
*
(
*
fct
)(
qp1
[
iq
],
qp2
[
iq
]);
value
+=
tmp
*
dualElInfo
.
smallElInfo
->
getDet
();
cont
=
dualTraverse
.
traverseNext
(
dualElInfo
);
}
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
double
localValue
=
value
;
MPI
::
COMM_WORLD
.
Allreduce
(
&
localValue
,
&
value
,
1
,
MPI_DOUBLE
,
MPI_SUM
);
#endif
return
value
;
}
double
integrate
(
const
DOFVector
<
double
>
&
vec
,
AbstractFunction
<
double
,
WorldVector
<
double
>
>
*
fct
)
{
FUNCNAME
(
"integrate()"
);
TEST_EXIT
(
fct
)(
"No function defined!
\n
"
);
const
FiniteElemSpace
*
feSpace
=
vec
.
getFeSpace
();
Mesh
*
mesh
=
feSpace
->
getMesh
();
int
deg
=
std
::
max
(
fct
->
getDegree
(),
2
*
feSpace
->
getBasisFcts
()
->
getDegree
());
Quadrature
*
quad
=
Quadrature
::
provideQuadrature
(
mesh
->
getDim
(),
deg
);
FastQuadrature
*
fastQuad
=
FastQuadrature
::
provideFastQuadrature
(
feSpace
->
getBasisFcts
(),
*
quad
,
INIT_PHI
);
mtl
::
dense_vector
<
double
>
qp
(
fastQuad
->
getNumPoints
());
WorldVector
<
double
>
coords
;
double
value
=
0.0
;
Flag
traverseFlag
=
Mesh
::
CALL_LEAF_EL
|
Mesh
::
FILL_COORDS
|
Mesh
::
FILL_DET
;
TraverseStack
stack
;
ElInfo
*
elInfo
=
stack
.
traverseFirst
(
mesh
,
-
1
,
traverseFlag
);
while
(
elInfo
)
{
vec
.
getVecAtQPs
(
elInfo
,
quad
,
fastQuad
,
qp
);
double
tmp
=
0.0
;
for
(
int
iq
=
0
;
iq
<
fastQuad
->
getNumPoints
();
iq
++
)
{
elInfo
->
coordToWorld
(
fastQuad
->
getLambda
(
iq
),
coords
);
tmp
+=
fastQuad
->
getWeight
(
iq
)
*
(
qp
[
iq
]
*
(
*
fct
)(
coords
));
}
value
+=
tmp
*
elInfo
->
getDet
();
elInfo
=
stack
.
traverseNext
(
elInfo
);
}
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
double
localValue
=
value
;
MPI
::
COMM_WORLD
.
Allreduce
(
&
localValue
,
&
value
,
1
,
MPI_DOUBLE
,
MPI_SUM
);
#endif
return
value
;
}
double
integrate
(
const
FiniteElemSpace
*
feSpace
,
AbstractFunction
<
double
,
WorldVector
<
double
>
>
*
fct
)
{
FUNCNAME
(
"integrate()"
);
TEST_EXIT
(
fct
)(
"No function defined!
\n
"
);
int
deg
=
std
::
max
(
fct
->
getDegree
(),
feSpace
->
getBasisFcts
()
->
getDegree
());
Quadrature
*
quad
=
Quadrature
::
provideQuadrature
(
feSpace
->
getMesh
()
->
getDim
(),
deg
);
FastQuadrature
*
fastQuad
=
FastQuadrature
::
provideFastQuadrature
(
feSpace
->
getBasisFcts
(),
*
quad
,
INIT_PHI
);
WorldVector
<
double
>
coords
;
double
value
=
0.0
;
Flag
traverseFlag
=
Mesh
::
CALL_LEAF_EL
|
Mesh
::
FILL_COORDS
|
Mesh
::
FILL_DET
;
TraverseStack
stack
;
ElInfo
*
elInfo
=
stack
.
traverseFirst
(
feSpace
->
getMesh
(),
-
1
,
traverseFlag
);
while
(
elInfo
)
{
double
tmp
=
0.0
;
for
(
int
iq
=
0
;
iq
<
fastQuad
->
getNumPoints
();
iq
++
)
{
elInfo
->
coordToWorld
(
fastQuad
->
getLambda
(
iq
),
coords
);
tmp
+=
fastQuad
->
getWeight
(
iq
)
*
(
*
fct
)(
coords
);
}
value
+=
tmp
*
elInfo
->
getDet
();
elInfo
=
stack
.
traverseNext
(
elInfo
);
}
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
double
localValue
=
value
;
MPI
::
COMM_WORLD
.
Allreduce
(
&
localValue
,
&
value
,
1
,
MPI_DOUBLE
,
MPI_SUM
);
#endif
return
value
;
}
double
integrateGeneral
(
const
std
::
vector
<
DOFVector
<
double
>*>
&
vecs
,
AbstractFunction
<
double
,
std
::
vector
<
double
>
>
*
fct
)
{
...
...
AMDiS/src/DOFVector.h
View file @
30debe8f
...
...
@@ -872,46 +872,89 @@ namespace AMDiS {
template
<
typename
T
>
std
::
vector
<
DOFVector
<
double
>*>
*
transform
(
DOFVector
<
typename
GradientType
<
T
>::
type
>
*
vec
,
std
::
vector
<
DOFVector
<
double
>*>
*
res
);
/// Computes the integral: \f$ \int f(\vec{x})\,\text{d}\vec{x}\f$
template
<
typename
TOut
>
TOut
integrate_Coords
(
const
FiniteElemSpace
*
feSpace
,
AbstractFunction
<
TOut
,
WorldVector
<
double
>
>
*
fct
);
template
<
typename
TOut
>
TOut
integrate
(
const
FiniteElemSpace
*
feSpace
,
AbstractFunction
<
TOut
,
WorldVector
<
double
>
>
*
fct
)
{
return
integrate_Coords
(
feSpace
,
fct
);
}
/// Computes the integral: \f$ \int f(v(\vec{x}))\,\text{d}\vec{x}\f$
template
<
typename
TOut
,
typename
T
>
TOut
integrate_Vec
(
const
DOFVector
<
T
>
&
vec
,
AbstractFunction
<
TOut
,
T
>
*
fct
);
template
<
typename
TOut
,
typename
T
>
TOut
integrate
(
const
DOFVector
<
T
>
&
vec
,
AbstractFunction
<
TOut
,
T
>
*
fct
)
{
return
integrate_Vec
(
vec
,
fct
);
}
/** \brief
* Computes the integral of a function that includes two different DOFVectors. This
* function works also for the case that the DOFVectors are defined on two different
* Computes the integral of a function that includes two different DOFVectors:
* \f$ \int f(v(\vec{x}), w(\vec{x}))\,\text{d}\vec{x}\f$
* This function works also for the case that the DOFVectors are defined on two different
* meshes.
*/
double
integrate
(
const
DOFVector
<
double
>
&
vec1
,
const
DOFVector
<
double
>
&
vec2
,
BinaryAbstractFunction
<
double
,
double
,
double
>
*
fct
);
template
<
typename
TOut
,
typename
T1
,
typename
T2
>
TOut
integrate_Vec2
(
const
DOFVector
<
T1
>
&
vec1
,
const
DOFVector
<
T2
>
&
vec2
,
BinaryAbstractFunction
<
TOut
,
T1
,
T2
>
*
fct
);
template
<
typename
TOut
,
typename
T1
,
typename
T2
>
TOut
integrate
(
const
DOFVector
<
T1
>
&
vec1
,
const
DOFVector
<
T2
>
&
vec2
,
BinaryAbstractFunction
<
TOut
,
T1
,
T2
>
*
fct
)
{
return
integrate_Vec2
(
vec1
,
vec2
,
fct
);
}
/// Computes the integral of a function with two DOFVectors defined on the same mesh.
double
intSingleMesh
(
const
DOFVector
<
double
>
&
vec1
,
const
DOFVector
<
double
>
&
vec2
,
BinaryAbstractFunction
<
double
,
double
,
double
>
*
fct
);
template
<
typename
TOut
,
typename
T1
,
typename
T2
>
TOut
int_Vec2_SingleMesh
(
const
DOFVector
<
T1
>
&
vec1
,
const
DOFVector
<
T2
>
&
vec2
,
BinaryAbstractFunction
<
TOut
,
T1
,
T2
>
*
fct
);
/// Computes the integral of a function with two DOFVectors defined on different meshes.
double
intDualMesh
(
const
DOFVector
<
double
>
&
vec1
,
const
DOFVector
<
double
>
&
vec2
,
BinaryAbstractFunction
<
double
,
double
,
double
>
*
fct
);
double
integrate
(
const
DOFVector
<
double
>
&
vec
,
AbstractFunction
<
double
,
WorldVector
<
double
>
>
*
fct
);
double
integrate
(
const
FiniteElemSpace
*
feSpace
,
AbstractFunction
<
double
,
WorldVector
<
double
>
>
*
fct
);
template
<
typename
TOut
,
typename
T1
,
typename
T2
>
TOut
int_Vec2_DualMesh
(
const
DOFVector
<
T1
>
&
vec1
,
const
DOFVector
<
T2
>
&
vec2
,
BinaryAbstractFunction
<
TOut
,
T1
,
T2
>
*
fct
);
/// Computes the integral: \f$ \int v(\vec{x}) f(\vec{x})\,\text{d}\vec{x}\f$
template
<
typename
T1
,
typename
T2
>
typename
ProductType
<
T1
,
T2
>::
type
integrate_VecTimesCoords
(
const
DOFVector
<
T1
>
&
vec
,
AbstractFunction
<
T2
,
WorldVector
<
double
>
>
*
fct
);
/// Computes the integral: \f$ \int f(v(\vec{x}), \vec{x})\,\text{d}\vec{x}\f$
template
<
typename
TOut
,
typename
T
>
TOut
integrate_VecAndCoords
(
const
DOFVector
<
T
>
&
vec
,
BinaryAbstractFunction
<
TOut
,
T
,
WorldVector
<
double
>
>
*
fct
);
template
<
typename
TOut
,
typename
T
>
TOut
integrate_VecAndCoords
(
const
DOFVector
<
T
>
&
vec
,
BinaryAbstractFunction
<
TOut
,
T
,
WorldVector
<
double
>
>
*
fct
)
{
return
integrate
(
vec
,
fct
);
}
/// Computes the integral: \f$ \int f(\{v_i(\vec{x})\}_i)\,\text{d}\vec{x}\f$
double
integrateGeneral
(
const
std
::
vector
<
DOFVector
<
double
>*>
&
vecs
,
AbstractFunction
<
double
,
std
::
vector
<
double
>
>
*
fct
);
/// Computes the integral: \f$ \int f(\{v_i(\vec{x})\}_i,\{\nabla w_i(\vec{x})\}_i)\,\text{d}\vec{x}\f$
double
integrateGeneralGradient
(
const
std
::
vector
<
DOFVector
<
double
>*>
&
vecs
,
const
std
::
vector
<
DOFVector
<
double
>*>
&
grds
,
BinaryAbstractFunction
<
double
,
std
::
vector
<
double
>
,
std
::
vector
<
WorldVector
<
double
>
>
>
*
fct
);
template
<
typename
T
>
T
integrate_VecAndCoords
(
const
DOFVector
<
double
>
&
vec
,
BinaryAbstractFunction
<
T
,
double
,
WorldVector
<
double
>
>
*
fct
);
template
<
typename
T1
,
typename
T2
>
T2
integrate_Vec
(
const
DOFVector
<
T1
>
&
vec
,
AbstractFunction
<
T2
,
T1
>
*
fct
);
}
#include
"DOFVector.hh"
...
...
AMDiS/src/DOFVector.hh
View file @
30debe8f
...
...
@@ -34,6 +34,11 @@
#include
"Operator.h"
#include
"Initfile.h"
#include
"Traverse.h"
#include
"DualTraverse.h"
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
#include
"MpiHelper.h"
#endif
// Defining the interface for MTL4
namespace
mtl
{
...
...
@@ -516,7 +521,7 @@ namespace AMDiS {
FastQuadrature
*
quadFast
=
FastQuadrature
::
provideFastQuadrature
(
this
->
feSpace
->
getBasisFcts
(),
*
q
,
INIT_PHI
);
double
result
=
0.0
;
T
result
;
nullify
(
result
)
;
int
nPoints
=
quadFast
->
getNumPoints
();
mtl
::
dense_vector
<
T
>
uh_vec
(
nPoints
);
TraverseStack
stack
;
...
...
@@ -524,7 +529,7 @@ namespace AMDiS {
Mesh
::
CALL_LEAF_EL
|
Mesh
::
FILL_COORDS
|
Mesh
::
FILL_DET
);
while
(
elInfo
)
{
double
det
=
elInfo
->
getDet
();
double
normT
=
0.0
;
T
normT
;
nullify
(
normT
)
;
this
->
getVecAtQPs
(
elInfo
,
NULL
,
quadFast
,
uh_vec
);
for
(
int
iq
=
0
;
iq
<
nPoints
;
iq
++
)
normT
+=
quadFast
->
getWeight
(
iq
)
*
(
uh_vec
[
iq
]);
...
...
@@ -534,18 +539,58 @@ namespace AMDiS {
}
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
double
localResult
=
result
;
MPI
::
COMM_WORLD
.
Allreduce
(
&
localResult
,
&
result
,
1
,
MPI_DOUBLE
,
MPI_SUM
);
mpi
::
globalAdd
(
result
);
#endif
return
result
;
}
template
<
typename
T
>
T
integrate_VecAndCoords
(
const
DOFVector
<
double
>
&
vec
,
BinaryAbstractFunction
<
T
,
double
,
WorldVector
<
double
>
>
*
fct
)
template
<
typename
TOut
>
TOut
integrate_Coords
(
const
FiniteElemSpace
*
feSpace
,
AbstractFunction
<
TOut
,
WorldVector
<
double
>
>
*
fct
)
{
FUNCNAME
(
"integrate_VecAndCoords()"
);
FUNCNAME
(
"integrate_Coords()"
);
TEST_EXIT
(
fct
)(
"No function defined!
\n
"
);
int
deg
=
std
::
max
(
fct
->
getDegree
(),
feSpace
->
getBasisFcts
()
->
getDegree
());
Quadrature
*
quad
=
Quadrature
::
provideQuadrature
(
feSpace
->
getMesh
()
->
getDim
(),
deg
);
FastQuadrature
*
fastQuad
=
FastQuadrature
::
provideFastQuadrature
(
feSpace
->
getBasisFcts
(),
*
quad
,
INIT_PHI
);
WorldVector
<
double
>
coords
;
TOut
value
;
nullify
(
value
);
Flag
traverseFlag
=
Mesh
::
CALL_LEAF_EL
|
Mesh
::
FILL_COORDS
|
Mesh
::
FILL_DET
;
TraverseStack
stack
;
ElInfo
*
elInfo
=
stack
.
traverseFirst
(
feSpace
->
getMesh
(),
-
1
,
traverseFlag
);
while
(
elInfo
)
{
TOut
tmp
;
nullify
(
tmp
);
for
(
int
iq
=
0
;
iq
<
fastQuad
->
getNumPoints
();
iq
++
)
{
elInfo
->
coordToWorld
(
fastQuad
->
getLambda
(
iq
),
coords
);
tmp
+=
fastQuad
->
getWeight
(
iq
)
*
(
*
fct
)(
coords
);
}
value
+=
tmp
*
elInfo
->
getDet
();
elInfo
=
stack
.
traverseNext
(
elInfo
);
}
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
mpi
::
globalAdd
(
value
);
#endif
return
value
;
}
template
<
typename
TOut
,
typename
T
>
TOut
integrate_Vec
(
const
DOFVector
<
T
>
&
vec
,
AbstractFunction
<
TOut
,
T
>
*
fct
)
{
FUNCNAME
(
"integrate_Vec()"
);
TEST_EXIT
(
fct
)(
"No function defined!
\n
"
);
...
...
@@ -556,39 +601,184 @@ namespace AMDiS {
FastQuadrature
*
fastQuad
=
FastQuadrature
::
provideFastQuadrature
(
vec
.
getFeSpace
()
->
getBasisFcts
(),
*
quad
,
INIT_PHI
);
mtl
::
dense_vector
<
double
>
qp
(
fastQuad
->
getNumPoints
());
WorldVector
<
double
>
coordsAtQP
;
mtl
::
dense_vector
<
T
>
qp
(
fastQuad
->
getNumPoints
());
TOut
value
;
nullify
(
value
);
T
value
;
nullify
(
value
);
Flag
traverseFlag
=
Mesh
::
CALL_LEAF_EL
|
Mesh
::
FILL_COORDS
|
Mesh
::
FILL_DET
;
TraverseStack
stack
;
ElInfo
*
elInfo
=
stack
.
traverseFirst
(
vec
.
getFeSpace
()
->
getMesh
(),
-
1
,
traverseFlag
);
while
(
elInfo
)
{
vec
.
getVecAtQPs
(
elInfo
,
quad
,
fastQuad
,
qp
);
T
tmp
;
nullify
(
tmp
);
TOut
tmp
;
nullify
(
tmp
);
for
(
int
iq
=
0
;
iq
<
fastQuad
->
getNumPoints
();
iq
++
)
{
elInfo
->
coordToWorld
(
fastQuad
->
getLambda
(
iq
),
coordsAtQP
);
tmp
+=
fastQuad
->
getWeight
(
iq
)
*
(
*
fct
)(
qp
[
iq
],
coordsAtQP
);
tmp
+=
fastQuad
->
getWeight
(
iq
)
*
(
*
fct
)(
qp
[
iq
]);
}
value
+=
tmp
*
elInfo
->
getDet
();
elInfo
=
stack
.
traverseNext
(
elInfo
);
}
// #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
// double localValue = value;
// MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
// #endif
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
mpi
::
globalAdd
(
value
);
#endif
return
value
;
}
template
<
typename
TOut
,
typename
T1
,
typename
T2
>
TOut
integrate_Vec2
(
const
DOFVector
<
T1
>
&
vec1
,
const
DOFVector
<
T2
>
&
vec2
,
BinaryAbstractFunction
<
TOut
,
T1
,
T2
>
*
fct
)
{
if
(
vec1
.
getFeSpace
()
->
getMesh
()
==
vec2
.
getFeSpace
()
->
getMesh
())
return
int_Vec2_SingleMesh
(
vec1
,
vec2
,
fct
);
else
return
int_Vec2_DualMesh
(
vec1
,
vec2
,
fct
);
}
template
<
typename
TOut
,
typename
T1
,
typename
T2
>
TOut
int_Vec2_SingleMesh
(
const
DOFVector
<
T1
>
&
vec1
,
const
DOFVector
<
T2
>
&
vec2
,
BinaryAbstractFunction
<
TOut
,
T1
,
T2
>
*
fct
)
{
FUNCNAME
(
"intDualmesh()"
);
TEST_EXIT
(
fct
)(
"No function defined!
\n
"
);
int
deg
=
std
::
max
(
fct
->
getDegree
(),
2
*
vec1
.
getFeSpace
()
->
getBasisFcts
()
->
getDegree
());
Quadrature
*
quad
=
Quadrature
::
provideQuadrature
(
vec1
.
getFeSpace
()
->
getMesh
()
->
getDim
(),
deg
);
FastQuadrature
*
fastQuad
=
FastQuadrature
::
provideFastQuadrature
(
vec1
.
getFeSpace
()
->
getBasisFcts
(),
*
quad
,
INIT_PHI
);
mtl
::
dense_vector
<
T1
>
qp1
(
fastQuad
->
getNumPoints
());
mtl
::
dense_vector
<
T2
>
qp2
(
fastQuad
->
getNumPoints
());
TOut
value
;
nullify
(
value
);
Flag
traverseFlag
=
Mesh
::
CALL_LEAF_EL
|
Mesh
::
FILL_COORDS
|
Mesh
::
FILL_DET
;
TraverseStack
stack
;
ElInfo
*
elInfo
=
stack
.
traverseFirst
(
vec1
.
getFeSpace
()
->
getMesh
(),
-
1
,
traverseFlag
);
while
(
elInfo
)
{
vec1
.
getVecAtQPs
(
elInfo
,
quad
,
fastQuad
,
qp1
);
vec2
.
getVecAtQPs
(
elInfo
,
quad
,
fastQuad
,
qp2
);
TOut
tmp
;
nullify
(
tmp
);
for
(
int
iq
=
0
;
iq
<
fastQuad
->
getNumPoints
();
iq
++
)
tmp
+=
fastQuad
->
getWeight
(
iq
)
*
(
*
fct
)(
qp1
[
iq
],
qp2
[
iq
]);
value
+=
tmp
*
elInfo
->
getDet
();
elInfo
=
stack
.
traverseNext
(
elInfo
);
}
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
mpi
::
globalAdd
(
value
);
#endif
return
value
;
}
template
<
typename
TOut
,
typename
T1
,
typename
T2
>
TOut
int_Vec2_DualMesh
(
const
DOFVector
<
T1
>
&
vec1
,
const
DOFVector
<
T2
>
&
vec2
,
BinaryAbstractFunction
<
TOut
,
T1
,
T2
>
*
fct
)
{
FUNCNAME
(
"intDualmesh()"
);
TEST_EXIT
(
fct
)(
"No function defined!
\n
"
);
int
deg
=
std
::
max
(
fct
->
getDegree
(),
2
*
vec1
.
getFeSpace
()
->
getBasisFcts
()
->
getDegree
());
Quadrature
*
quad
=
Quadrature
::
provideQuadrature
(
vec1
.
getFeSpace
()
->
getMesh
()
->
getDim
(),
deg
);
FastQuadrature
*
fastQuad
=
FastQuadrature
::
provideFastQuadrature
(
vec1
.
getFeSpace
()
->
getBasisFcts
(),
*
quad
,
INIT_PHI
);
mtl
::
dense_vector
<
T1
>
qp1
(
fastQuad
->
getNumPoints
());
mtl
::
dense_vector
<
T2
>
qp2
(
fastQuad
->
getNumPoints
());
TOut
value
;
nullify
(
value
);
Flag
traverseFlag
=
Mesh
::
CALL_LEAF_EL
|
Mesh
::
FILL_COORDS
|
Mesh
::
FILL_DET
;
DualTraverse
dualTraverse
;
DualElInfo
dualElInfo
;
bool
cont
=
dualTraverse
.
traverseFirst
(
vec1
.
getFeSpace
()
->
getMesh
(),
vec2
.
getFeSpace
()
->
getMesh
(),
-
1
,
-
1
,
traverseFlag
,
traverseFlag
,
dualElInfo
);
while
(
cont
)
{
vec1
.
getVecAtQPs
(
dualElInfo
.
smallElInfo
,
dualElInfo
.
largeElInfo
,
quad
,
NULL
,
qp1
);
vec2
.
getVecAtQPs
(
dualElInfo
.
smallElInfo
,
dualElInfo
.
largeElInfo
,
quad
,
NULL
,
qp2
);
TOut
tmp
;
nullify
(
tmp
);
for
(
int
iq
=
0
;
iq
<
quad
->
getNumPoints
();
iq
++
)
tmp
+=
quad
->
getWeight
(
iq
)
*
(
*
fct
)(
qp1
[
iq
],
qp2
[
iq
]);
value
+=
tmp
*
dualElInfo
.
smallElInfo
->
getDet
();
cont
=
dualTraverse
.
traverseNext
(
dualElInfo
);
}
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
mpi
::
globalAdd
(
value
);
#endif
return
value
;
}
template
<
typename
T1
,
typename
T2
>
T2
integrate_Vec
(
const
DOFVector
<
T1
>
&
vec
,
AbstractFunction
<
T2
,
T1
>
*
fct
)
typename
ProductType
<
T1
,
T2
>::
type
integrate_VecTimesCoords
(
const
DOFVector
<
T1
>
&
vec
,
AbstractFunction
<
T2
,
WorldVector
<
double
>
>
*
fct
)
{
FUNCNAME
(
"integrate_VecTimesCoords()"
);
TEST_EXIT
(
fct
)(
"No function defined!
\n
"
);
typedef
typename
ProductType
<
T1
,
T2
>::
type
TOut
;
const
FiniteElemSpace
*
feSpace
=
vec
.
getFeSpace
();
Mesh
*
mesh
=
feSpace
->
getMesh
();
int
deg
=
std
::
max
(
fct
->
getDegree
(),
2
*
feSpace
->
getBasisFcts
()
->
getDegree
());
Quadrature
*
quad
=
Quadrature
::
provideQuadrature
(
mesh
->
getDim
(),
deg
);
FastQuadrature
*
fastQuad
=
FastQuadrature
::
provideFastQuadrature
(
feSpace
->
getBasisFcts
(),
*
quad
,
INIT_PHI
);
mtl
::
dense_vector
<
T1
>
qp
(
fastQuad
->
getNumPoints
());
WorldVector
<
double
>
coords
;
TOut
value
;
nullify
(
value
);
Flag
traverseFlag
=
Mesh
::
CALL_LEAF_EL
|
Mesh
::
FILL_COORDS
|
Mesh
::
FILL_DET
;
TraverseStack
stack
;
ElInfo
*
elInfo
=
stack
.
traverseFirst
(
mesh
,
-
1
,
traverseFlag
);
while
(
elInfo
)
{
vec
.
getVecAtQPs
(
elInfo
,
quad
,
fastQuad
,
qp
);
TOut
tmp
;
nullify
(
tmp
);
for
(
int
iq
=
0
;
iq
<
fastQuad
->
getNumPoints
();
iq
++
)
{
elInfo
->
coordToWorld
(
fastQuad
->
getLambda
(
iq
),
coords
);
tmp
+=
fastQuad
->
getWeight
(
iq
)
*
(
qp
[
iq
]
*
(
*
fct
)(
coords
));
}
value
+=