Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in
Toggle navigation
Menu
Open sidebar
Mpp
MLUQ
Commits
71788f0f
Commit
71788f0f
authored
Aug 10, 2020
by
niklas.baumgarten
Browse files
new strucutre (major refactoring like in tutorial)
parent
63c4337a
Pipeline
#103787
failed with stages
in 12 minutes and 4 seconds
Changes
62
Pipelines
1
Hide whitespace changes
Inline
Side-by-side
mlmc/src/CMakeLists.txt
View file @
71788f0f
set
(
MLMC_SRC
Main.cpp
utils/Utils.cpp
MainProgram.cpp
MultilevelMonteCarlo.cpp
mc/MonteCarlo.cpp
mc/MonteCarloElliptic.cpp
mc/MonteCarloTransport.cpp
assemble/EllipticAssemble.C
assemble/MixedEllipticAssemble.C
assemble/HybridEllipticAssemble.C
assemble/DGEllipticAssemble.C
stochastics/RNManager.C
stochastics/CirculantEmbedding.C
problem/StochasticEllipticProblem.hpp
stochastics/SampleGenerator.h
stochastics/HybridFluxGenerator.C
stochastics/HybridFluxGenerator.h
mc/MonteCarlo.cpp
utils/MultilevelPlotter.cpp
utils/MultilevelPlotter.hpp
)
include
(
assemble/CMakeLists.txt
)
include
(
main/CMakeLists.txt
)
include
(
montecarlo/CMakeLists.txt
)
include
(
problem/CMakeLists.txt
)
include
(
stochastics/CMakeLists.txt
)
add_library
(
MLMC STATIC
${
MLMC_SRC
}
)
add_library
(
MLMC STATIC
${
assemble
}
${
main
}
${
montecarlo
}
${
problem
}
${
stochastics
}
)
mlmc/src/Main.cpp
View file @
71788f0f
#include
"m++.hpp"
#include
"MainProgram.hpp"
#include
"
main/
MainProgram.hpp"
#ifdef COMPLEX
#error undef COMPLEX in src/Compiler.h
...
...
mlmc/src/MonteCarloLogger.C
deleted
100644 → 0
View file @
63c4337a
#include
"MonteCarloLogger.h"
IndentedLogger
*
IndentedLogger
::
instance
=
nullptr
;
mlmc/src/MonteCarloLogger.h
deleted
100644 → 0
View file @
63c4337a
#ifndef MLMC__MULTILEVELMONTECARLO__HPP
#define MLMC__MULTILEVELMONTECARLO__HPP
class
IndentedLogger
{
static
IndentedLogger
*
instance
;
string
indent
;
string
innerIndent
;
const
string
defaultIndent
=
" "
;
vector
<
Date
>
timestamps
;
IndentedLogger
()
{
indent
=
""
;
}
public:
static
IndentedLogger
*
GetInstance
()
{
if
(
!
instance
)
instance
=
new
IndentedLogger
;
return
instance
;
}
int
DecreaseIndent
()
{
indent
=
indent
.
erase
(
0
,
defaultIndent
.
length
());
}
void
IncreaseIndent
()
{
indent
=
indent
.
append
(
defaultIndent
);
}
string
GetIndent
()
{
return
indent
;
}
void
StartMethod
(
const
string
&
msg
,
int
verbose
)
{
Date
start
;
timestamps
.
push_back
(
start
);
vout
(
1
)
<<
indent
<<
"<"
<<
msg
<<
">"
<<
endl
;
IncreaseIndent
();
}
void
EndMethod
(
int
verbose
)
{
DecreaseIndent
();
vout
(
1
)
<<
indent
<<
"<End after "
<<
Date
()
-
timestamps
.
back
()
<<
">"
<<
endl
;
timestamps
.
erase
(
timestamps
.
end
());
}
void
InnerMsg
(
const
string
&
msg
,
int
verbose
)
{
vout
(
1
)
<<
indent
<<
msg
<<
endl
;
}
void
LogMsgFlush
(
const
string
&
msg
,
int
verbose
)
{
vout
(
1
)
<<
"
\r
"
<<
indent
<<
msg
<<
flush
;
}
};
#endif //MLMC__MULTILEVELMONTECARLO__HPP
mlmc/src/assemble/AssembleHeader.h
deleted
100644 → 0
View file @
63c4337a
#ifndef MLMC_ASSEMBLEHEADER_H
#define MLMC_ASSEMBLEHEADER_H
#include
"EllipticAssemble.h"
#include
"MixedEllipticAssemble.h"
#include
"HybridEllipticAssemble.h"
#include
"DGEllipticAssemble.h"
#include
"DGTransportAssemble.h"
#endif //MLMC_ASSEMBLEHEADER_H
mlmc/src/assemble/CMakeLists.txt
0 → 100644
View file @
71788f0f
set
(
assemble
assemble/EllipticAssemble.cpp
assemble/MixedEllipticAssemble.cpp
assemble/DGEllipticAssemble.cpp
assemble/HybridEllipticAssemble.cpp
assemble/DGTransportAssemble.cpp
assemble/DGReactionAssemble.cpp
assemble/PGReactionAssemble.cpp
)
\ No newline at end of file
mlmc/src/assemble/DGEllipticAssemble.
C
→
mlmc/src/assemble/DGEllipticAssemble.
cpp
View file @
71788f0f
#include
"DGEllipticAssemble.h"
Element
*
DGEllipticAssemble
::
getElement
(
const
cell
&
c
,
const
Vector
&
u
)
const
{
return
new
DGElement
(
*
disc
,
u
,
c
);
}
FaceElement
*
DGEllipticAssemble
::
getFaceElement
(
const
cell
&
c
,
const
Vector
&
u
,
int
face
)
const
{
return
new
DGFaceElement
(
*
disc
,
u
,
c
,
face
);
}
Scalar
DGEllipticAssemble
::
getValue
(
const
cell
&
c
,
const
Vector
&
u
,
Element
*
elem
,
int
q
)
const
{
auto
dgEllipticElem
=
dynamic_cast
<
DGElement
*>
(
elem
);
return
dgEllipticElem
->
Value
(
q
,
u
);
}
Scalar
DGEllipticAssemble
::
getFaceValue
(
const
cell
&
c
,
const
Vector
&
u
,
FaceElement
*
faceElem
,
int
q
,
int
face
)
const
{
auto
dgEllipticFaceElem
=
dynamic_cast
<
DGFaceElement
*>
(
faceElem
);
return
dgEllipticFaceElem
->
Value
(
q
,
u
);
}
VectorField
DGEllipticAssemble
::
getFlux
(
const
cell
&
c
,
const
Vector
&
u
,
Element
*
elem
,
int
q
)
const
{
auto
dgEllipticElem
=
dynamic_cast
<
DGElement
*>
(
elem
);
Tensor
K
=
problem
->
Permeability
(
c
);
return
K
*
dgEllipticElem
->
Derivative
(
q
,
u
);
}
VectorField
DGEllipticAssemble
::
getFaceFlux
(
const
cell
&
c
,
const
Vector
&
u
,
FaceElement
*
faceElem
,
int
q
,
int
face
)
const
{
auto
dgEllipticFaceElem
=
dynamic_cast
<
DGFaceElement
*>
(
faceElem
);
Tensor
K
=
problem
->
Permeability
(
c
);
return
K
*
dgEllipticFaceElem
->
Derivative
(
q
,
u
);
}
VectorField
DGEllipticAssemble
::
getGradient
(
const
cell
&
c
,
const
Vector
&
u
,
Element
*
elem
,
int
q
)
const
{
auto
dgEllipticElem
=
dynamic_cast
<
DGElement
*>
(
elem
);
return
dgEllipticElem
->
Derivative
(
q
,
u
);
}
#include
"DGEllipticAssemble.hpp"
void
DGEllipticAssemble
::
Dirichlet
(
const
cell
&
c
,
Vector
&
u
)
const
{}
...
...
@@ -218,6 +176,19 @@ void DGEllipticAssemble::Jacobi(const cell &c, const Vector &u, Matrix &A) const
}
}
void
DGEllipticAssemble
::
SetExactSolution
(
Vector
&
u
)
const
{
u
=
0
;
Point
z
[
MaxNodalPoints
];
for
(
cell
c
=
u
.
cells
();
c
!=
u
.
cells_end
();
++
c
)
{
row
r
=
u
.
find_row
(
c
());
DGElement
Elem
(
*
disc
,
u
,
c
);
Elem
.
NodalPoints
(
z
);
for
(
int
j
=
0
;
j
<
Elem
.
dimension
();
++
j
)
u
(
r
,
j
)
=
problem
->
Solution
(
z
[
j
]);
}
Accumulate
(
u
);
}
void
DGEllipticAssemble
::
AssembleTransfer
(
TransferMatrix
&
TM
)
const
{
TM
=
0
;
const
matrixgraph
&
cg
=
TM
.
CoarseMatrixGraph
();
...
...
@@ -228,22 +199,20 @@ void DGEllipticAssemble::AssembleTransfer(TransferMatrix &TM) const {
TM
(
c
(),
c
.
Child
(
k
))[
0
]
=
1
;
}
//void plotU(const Vector &u, const string &name) override;
//
//void DGEllipticAssemble::plotU(const Vector &u, const string &name) {
// Vector plot_u(u);
// for (cell c = u.cells(); c != u.cells_end(); ++c) {
// DGElement elem(*disc, u, c);
// double U = 0;
// double area = 0;
// for (int q = 0; q < elem.nQ(); ++q) {
// double w = elem.QWeight(q);
// U += w * elem.Value(q, u);
// area += w;
// }
// U *= (1 / area);
// plot_u(c(), 0) = U;
// }
// plot->celldata(plot_u);
// plot->vtk_celldata("u");
//}
void
DGEllipticAssemble
::
plotU
(
const
Vector
&
u
,
Meshes
&
meshes
)
{
Vector
plot_u
(
u
);
for
(
cell
c
=
u
.
cells
();
c
!=
u
.
cells_end
();
++
c
)
{
DGElement
elem
(
*
disc
,
u
,
c
);
double
U
=
0
;
double
area
=
0
;
for
(
int
q
=
0
;
q
<
elem
.
nQ
();
++
q
)
{
double
w
=
elem
.
QWeight
(
q
);
U
+=
w
*
elem
.
Value
(
q
,
u
);
area
+=
w
;
}
U
*=
(
1
/
area
);
plot_u
(
c
(),
0
)
=
U
;
}
// plot->celldata(plot_u);
// plot->vtk_celldata("u_ex");
}
mlmc/src/assemble/DGEllipticAssemble.h
→
mlmc/src/assemble/DGEllipticAssemble.h
pp
View file @
71788f0f
#ifndef _DGLAPLACE_H_
#define _DGLAPLACE_H_
#include
"EllipticAssemble.h"
#include
"EllipticAssemble.h
pp
"
#include
"discretization/DGDiscretization.hpp"
#include
"discretization/DGUtils.h"
class
DGEllipticAssemble
:
public
EllipticAssemble
,
public
DGUtils
{
public:
std
::
shared_ptr
<
DGDiscretization
>
disc
;
DGDiscretization
*
disc
;
double
penalty
=
1.0
;
int
sign
=
1
;
DGEllipticAssemble
(
std
::
shared_ptr
<
IDiscretizationT
<>>
disc
,
std
::
shared_ptr
<
StochasticEllipticProblem
>
problem
)
:
EllipticAssemble
(
move
(
disc
),
move
(
problem
))
{
this
->
disc
=
std
::
dynamic_pointer_cast
<
DGDiscretization
>
(
disc
);
DGEllipticAssemble
(
IDiscretization
*
disc
,
StochasticEllipticProblem
*
problem
)
:
EllipticAssemble
(
disc
,
problem
)
{
this
->
disc
=
dynamic_cast
<
DGDiscretization
*>
(
disc
);
config
.
get
(
"penalty"
,
penalty
);
config
.
get
(
"sign"
,
sign
);
}
const
char
*
Name
()
const
override
{
return
"Discontinuous-Galerkin"
;
}
Element
*
getElement
(
const
cell
&
c
,
const
Vector
&
u
)
const
override
;
Element
*
createElement
(
const
cell
&
c
,
const
Vector
&
u
)
const
override
{
return
new
DGElement
(
*
disc
,
u
,
c
);
}
FaceElement
*
getFaceElement
(
const
cell
&
c
,
const
Vector
&
u
,
int
face
)
const
override
;
FaceElement
*
createFaceElement
(
const
cell
&
c
,
const
Vector
&
u
,
int
face
)
const
override
{
return
new
DGFaceElement
(
*
disc
,
u
,
c
,
face
);
}
Scalar
getValue
(
const
cell
&
c
,
const
Vector
&
u
,
Element
*
elem
,
int
q
)
const
override
;
Scalar
getValue
(
const
cell
&
c
,
const
Vector
&
u
,
Element
*
elem
,
int
q
)
const
override
{
auto
dgEllipticElem
=
dynamic_cast
<
DGElement
*>
(
elem
);
return
dgEllipticElem
->
Value
(
q
,
u
);
}
Scalar
getFaceValue
(
const
cell
&
c
,
const
Vector
&
u
,
FaceElement
*
faceElem
,
int
q
,
int
face
)
const
override
;
FaceElement
*
faceElem
,
int
q
,
int
face
)
const
override
{
auto
dgEllipticFaceElem
=
dynamic_cast
<
DGFaceElement
*>
(
faceElem
);
return
dgEllipticFaceElem
->
Value
(
q
,
u
);
}
VectorField
getFlux
(
const
cell
&
c
,
const
Vector
&
u
,
Element
*
elem
,
int
q
)
const
override
;
Element
*
elem
,
int
q
)
const
override
{
auto
dgEllipticElem
=
dynamic_cast
<
DGElement
*>
(
elem
);
Tensor
K
=
problem
->
Permeability
(
c
);
return
K
*
dgEllipticElem
->
Derivative
(
q
,
u
);
}
VectorField
getFaceFlux
(
const
cell
&
c
,
const
Vector
&
u
,
FaceElement
*
faceElem
,
int
q
,
int
face
)
const
override
;
FaceElement
*
faceElem
,
int
q
,
int
face
)
const
override
{
auto
dgEllipticFaceElem
=
dynamic_cast
<
DGFaceElement
*>
(
faceElem
);
Tensor
K
=
problem
->
Permeability
(
c
);
return
K
*
dgEllipticFaceElem
->
Derivative
(
q
,
u
);
}
VectorField
getGradient
(
const
cell
&
c
,
const
Vector
&
u
,
Element
*
elem
,
int
q
)
const
override
;
Element
*
elem
,
int
q
)
const
override
{
auto
dgEllipticElem
=
dynamic_cast
<
DGElement
*>
(
elem
);
return
dgEllipticElem
->
Derivative
(
q
,
u
);
}
void
Dirichlet
(
const
cell
&
c
,
Vector
&
u
)
const
override
;
...
...
@@ -47,7 +70,13 @@ public:
void
Jacobi
(
const
cell
&
c
,
const
Vector
&
u
,
Matrix
&
A
)
const
override
;
void
SetExactSolution
(
Vector
&
u
)
const
override
;
void
AssembleTransfer
(
TransferMatrix
&
TM
)
const
override
;
void
plotU
(
const
Vector
&
u
,
Meshes
&
meshes
);
void
plotUex
(
const
Vector
&
u
,
Meshes
&
meshes
);
};
#endif
mlmc/src/assemble/DGReactionAssemble.cpp
0 → 100644
View file @
71788f0f
//
// Created by niklas on 10.08.20.
//
mlmc/src/assemble/DGReactionAssemble.h
deleted
100644 → 0
View file @
63c4337a
#ifndef _DG_REACTION_H_
#define _DG_REACTION_H_
#include
"dg/DGTAssemble.h"
class
DGReactionAssemble
:
public
DGTAssemble
{
public:
const
ScalarProblem
&
problem
;
double
penalty
=
1.0
;
int
sign
=
1
;
DGReactionAssemble
(
const
DGDiscretization
&
disc
,
const
ScalarProblem
&
problem
,
Plot
&
plot
)
:
DGTAssemble
(
disc
,
plot
),
problem
(
problem
)
{
ReadConfig
(
Settings
,
"penalty"
,
penalty
);
ReadConfig
(
Settings
,
"sign"
,
sign
);
}
const
char
*
Name
()
const
override
{
return
"DGReactionAssemble"
;
}
void
PrintInfo
(
const
Vector
&
u
)
override
{
Assemble
::
PrintInfo
(
u
);
Vout
(
1
)
<<
" Problem: "
<<
problem
.
Name
()
<<
endl
<<
" Discretization: "
<<
disc
.
DiscName
()
<<
endl
<<
endl
;
}
double
Residual
(
double
t
,
const
Vector
&
u
,
Vector
&
r
)
const
override
{
r
=
0
;
const
Vector
&
u_old
=
U_old
();
for
(
cell
c
=
u
.
cells
();
c
!=
u
.
cells_end
();
++
c
)
{
int
sd
=
c
.
Subdomain
();
RowBndValues
r_c
(
r
,
c
);
DGElement
elem
(
disc
,
u
,
c
);
double
kappa
=
0
;
for
(
int
q
=
0
;
q
<
elem
.
nQ
();
++
q
)
{
double
w
=
elem
.
QWeight
(
q
);
Point
z
=
elem
.
QPoint
(
q
);
Tensor
K
=
problem
.
Diffusion
(
sd
,
z
);
kappa
=
max
(
kappa
,
norm
(
K
));
Scalar
f
=
problem
.
Load
(
sd
,
z
);
Scalar
U
=
elem
.
Value
(
q
,
u
);
VectorField
DU
=
elem
.
Derivative
(
q
,
u
);
VectorField
K_DU
=
K
*
DU
;
Scalar
Uold
=
elem
.
Value
(
q
,
u_old
);
Scalar
R
=
problem
.
Reaction
(
z
,
U
);
VectorField
B
=
problem
.
cellConvection
(
sd
,
c
,
z
);
for
(
int
i
=
0
;
i
<
elem
.
dimension
();
++
i
)
{
Scalar
U_i
=
elem
.
Value
(
q
,
i
);
VectorField
DU_i
=
elem
.
Derivative
(
q
,
i
);
r
(
c
(),
i
)
+=
w
*
(
K_DU
*
DU_i
+
(
B
*
DU
)
*
U_i
+
(
1.0
/
dt_
)
*
(
U
-
Uold
)
*
U_i
-
(
R
+
f
)
*
U_i
);
}
}
for
(
int
f
=
0
;
f
<
c
.
Faces
();
++
f
)
{
Scalar
s
=
penalty
*
c
.
FaceArea
(
f
);
DGFaceElement
felem
(
disc
,
u
,
c
,
f
);
if
(
checkIfBND
(
r
,
c
,
f
))
{
for
(
int
q
=
0
;
q
<
felem
.
nQ
();
++
q
)
{
double
w
=
felem
.
QWeight
(
q
);
Point
z
=
felem
.
QPoint
(
q
);
Point
N
=
felem
.
QNormal
(
q
);
Scalar
BN
=
problem
.
faceConvection
(
sd
,
c
,
f
,
N
,
z
);
Tensor
K
=
problem
.
Diffusion
(
c
.
Subdomain
(),
z
);
Scalar
U_diff
=
felem
.
Value
(
q
,
u
)
-
problem
.
Concentration
(
t
,
z
);
VectorField
K_DU
=
K
*
felem
.
Derivative
(
q
,
u
);
if
(
BN
<
0
)
{
Scalar
F
=
K_DU
*
N
;
for
(
int
i
=
0
;
i
<
felem
.
dimension
();
++
i
)
{
Scalar
phi_i
=
felem
.
Value
(
q
,
i
);
Scalar
NDphi_i
=
(
K
*
felem
.
Derivative
(
q
,
i
))
*
N
;
r
(
c
(),
i
)
-=
w
*
(
F
*
phi_i
+
BN
*
U_diff
*
phi_i
+
U_diff
*
(
NDphi_i
-
s
*
phi_i
));
}
}
else
{
Scalar
F
=
problem
.
Flux
(
2
,
z
)
*
N
;
for
(
int
i
=
0
;
i
<
felem
.
dimension
();
++
i
)
{
Scalar
phi_i
=
felem
.
Value
(
q
,
i
);
r
(
c
(),
i
)
-=
w
*
F
*
phi_i
;
}
}
}
}
else
{
cell
cf
=
findNBCell
(
r
,
c
,
f
);
if
(
cf
()
<
c
())
continue
;
int
f1
=
findNBFaceID
(
r
,
c
.
Face
(
f
),
cf
);
DGFaceElement
felem_1
(
disc
,
r
,
cf
,
f1
);
for
(
int
q
=
0
;
q
<
felem
.
nQ
();
++
q
)
{
double
w
=
felem
.
QWeight
(
q
);
Point
z
=
felem
.
QPoint
(
q
);
Point
N
=
felem
.
QNormal
(
q
);
Scalar
BN
=
problem
.
faceConvection
(
sd
,
c
,
f
,
N
,
z
);
Tensor
K
=
problem
.
Diffusion
(
c
.
Subdomain
(),
z
);
Scalar
U
=
felem
.
Value
(
q
,
u
);
VectorField
K_DU
=
K
*
felem
.
Derivative
(
q
,
u
);
Point
Qf_c
=
felem
.
QPoint
(
q
);
int
q1
=
findQPointID
(
felem
,
felem_1
,
Qf_c
);
Tensor
K_1
=
problem
.
Diffusion
(
cf
.
Subdomain
(),
z
);
Scalar
U_1
=
felem_1
.
Value
(
q1
,
u
);
VectorField
K_DU_1
=
K_1
*
felem_1
.
Derivative
(
q1
,
u
);
for
(
int
i
=
0
;
i
<
felem
.
dimension
();
++
i
)
{
Scalar
phi_i
=
felem
.
Value
(
q
,
i
);
Scalar
NDphi_i
=
(
K
*
felem
.
Derivative
(
q
,
i
))
*
N
;
r
(
c
(),
i
)
+=
w
*
((
s
*
U
-
0.5
*
(
K_DU
*
N
))
*
phi_i
-
0.5
*
U
*
NDphi_i
);
r
(
c
(),
i
)
+=
w
*
((
-
s
*
U_1
-
0.5
*
(
K_DU_1
*
N
))
*
phi_i
+
0.5
*
U_1
*
NDphi_i
);
if
(
BN
<
0
)
{
r
(
c
(),
i
)
-=
w
*
BN
*
U
*
phi_i
;
r
(
c
(),
i
)
+=
w
*
BN
*
U_1
*
phi_i
;
}
}
for
(
int
j
=
0
;
j
<
felem_1
.
dimension
();
++
j
)
{
Scalar
phi_j
=
felem_1
.
Value
(
q1
,
j
);
Scalar
NDphi_j
=
(
K_1
*
felem_1
.
Derivative
(
q1
,
j
))
*
N
;
r
(
cf
(),
j
)
+=
w
*
((
0.5
*
K_DU_1
*
N
+
s
*
U_1
)
*
phi_j
+
0.5
*
U_1
*
NDphi_j
);
r
(
cf
(),
j
)
+=
w
*
((
0.5
*
K_DU
*
N
-
s
*
U
)
*
phi_j
-
0.5
*
U
*
NDphi_j
);
if
(
BN
>
0
)
{
r
(
cf
(),
j
)
+=
w
*
BN
*
U_1
*
phi_j
;
r
(
cf
(),
j
)
-=
w
*
BN
*
U
*
phi_j
;
}
}
}
}
}
}
Collect
(
r
);
return
r
.
norm
();
}
void
Jacobi
(
double
t
,
const
Vector
&
u
,
Matrix
&
A
)
const
override
{
A
=
0
;
for
(
cell
c
=
A
.
cells
();
c
!=
A
.
cells_end
();
++
c
)
{
int
sd
=
c
.
Subdomain
();
DGElement
elem
(
disc
,
A
,
c
);
DGRowEntries
A_c
(
A
,
c
,
c
);
double
kappa
=
0
;
for
(
int
q
=
0
;
q
<
elem
.
nQ
();
++
q
)
{
double
w
=
elem
.
QWeight
(
q
);
Point
z
=
elem
.
QPoint
(
q
);
Tensor
K
=
problem
.
Diffusion
(
sd
,
z
);
kappa
=
max
(
kappa
,
norm
(
K
));
Scalar
U
=
elem
.
Value
(
q
,
u
);
Scalar
DR
=
problem
.
D_Reaction
(
z
,
U
);
VectorField
B
=
problem
.
cellConvection
(
sd
,
c
,
z
);
for
(
int
i
=
0
;
i
<
elem
.
dimension
();
++
i
)
{
Scalar
U_i
=
elem
.
Value
(
q
,
i
);
VectorField
DU_i
=
elem
.
Derivative
(
q
,
i
);
VectorField
K_DU_i
=
K
*
DU_i
;
for
(
int
j
=
0
;
j
<
elem
.
dimension
();
++
j
)
{
Scalar
U_j
=
elem
.
Value
(
q
,
j
);
VectorField
DU_j
=
elem
.
Derivative
(
q
,
j
);
A_c
(
i
,
j
)
+=
w
*
(
K_DU_i
*
DU_j
+
(
B
*
DU_j
)
*
U_i
+
(
1.0
/
dt_
)
*
U_i
*
U_j
-
DR
*
U_i
*
U_j
);
}
}
}
BFParts
bnd
(
u
.
GetMesh
(),
c
);
for
(
int
f
=
0
;
f
<
c
.
Faces
();
++
f
)
{
Scalar
s
=
penalty
*
c
.
FaceArea
(
f
);
DGFaceElement
felem
(
disc
,
u
,
c
,
f
);
if
(
checkIfBND
(
u
,
c
,
f
))
{
for
(
int
q
=
0
;
q
<
felem
.
nQ
();
++
q
)
{
double
w
=
felem
.
QWeight
(
q
);
Point
z
=
felem
.
QPoint
(
q
);
Point
N
=
felem
.
QNormal
(
q
);
Scalar
BN
=
problem
.
faceConvection
(
sd
,
c
,
f
,
N
,
z
);
if
(
BN
<
0
)
{
Tensor
K
=
problem
.
Diffusion
(
c
.
Subdomain
(),
z
);
for
(
int
i
=
0
;
i
<
felem
.
dimension
();
++
i
)
{
Scalar
phi_i
=
felem
.
Value
(
q
,
i
);
Scalar
NDphi_i
=
(
K
*
felem
.
Derivative
(
q
,
i
))
*
N
;
for
(
int
j
=
0
;
j
<
felem
.
dimension
();
++
j
)
{
Scalar
phi_j
=
felem
.
Value
(
q
,
j
);
Scalar
NDphi_j
=
(
K
*
felem
.
Derivative
(
q
,
j
))
*
N
;
A_c
(
i
,
j
)
-=
w
*
(
NDphi_i
*
phi_j
+
BN
*
phi_j
*
phi_i
+
phi_i
*
(
NDphi_j
-
s
*
phi_j
));
}
}
}
}
}
else
{
cell
cf
=
findNBCell
(
u
,
c
,
f
);
if
(
cf
()
<
c
())
continue
;
DGRowEntries
A_cf
(
A
,
c
,
cf
);
DGRowEntries
A_fc
(
A
,
cf
,
c
);
DGRowEntries
A_ff
(
A
,
cf
,
cf
);
int
f1
=
findNBFaceID
(
u
,
c
.
Face
(
f
),
cf
);
DGFaceElement
felem_1
(
disc
,
u
,
cf
,
f1
);
for
(
int
q
=
0
;
q
<
felem
.
nQ
();
++
q
)
{
double
w
=
felem
.
QWeight
(
q
);
Point
z
=
felem
.
QPoint
(
q
);
Point
N
=
felem
.
QNormal
(
q
);
Scalar
BN
=
problem
.
faceConvection
(
sd
,
c
,
f
,
N
,
z
);
Tensor
K
=
problem
.
Diffusion
(
c
.
Subdomain
(),
z
);
Point
Qf_c
=
felem
.
QPoint
(
q
);
int
q1
=
findQPointID
(
felem
,
felem_1
,
Qf_c
);
Tensor
K_1
=
problem
.
Diffusion
(
cf
.
Subdomain
(),
z
);
for
(
int
i
=
0
;
i
<
felem
.
dimension
();
++
i
)
{
Scalar
phi_i
=
felem
.
Value
(
q
,
i
);
Scalar
NDphi_i
=
(
K
*
felem
.
Derivative
(
q
,
i
))
*
N
;
for
(
int
j
=
0
;
j
<
felem
.
dimension
();
++
j
)
{
Scalar
phi_j
=
felem
.
Value
(
q
,
j
);
Scalar
NDphi_j
=
(
K
*
felem
.
Derivative
(
q
,
j
))
*
N
;
A_c
(
i
,
j
)
+=
w
*
(
-
0.5
*
NDphi_i
*
phi_j
-
0.5
*
phi_i
*
NDphi_j
+
s
*
phi_i
*
phi_j
);
if
(
BN
<
0
)
A_c
(
i
,
j
)
-=
w
*
BN
*
phi_j
*
phi_i
;
}
for
(
int
j
=
0
;
j
<
felem_1
.
dimension
();
++
j
)
{
Scalar
phi_j
=
felem_1
.
Value
(
q1
,
j
);
Scalar
NDphi_j
=
(
K_1
*
felem_1
.
Derivative
(
q1
,
j
))
*
N
;
A_cf
(
i
,
j
)
+=
w
*
(
0.5
*
NDphi_i
*
phi_j
-
phi_i
*
0.5
*
NDphi_j
-
s
*
phi_i
*
phi_j
);
A_fc
(
j
,
i
)
+=
w
*
(
0.5
*
NDphi_i
*
phi_j
-
phi_i
*
0.5
*
NDphi_j
-
s
*
phi_i
*
phi_j
);
if
(
BN
<
0
)
A_cf
(
i
,
j
)
+=
w
*
BN
*
phi_j
*
phi_i
;
else
A_fc
(
j
,
i
)
-=
w
*
BN
*
phi_j
*
phi_i
;