Skip to content
GitLab
Menu
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in
Toggle navigation
Menu
Open sidebar
Tilman Diego Aleman
reproduction
Commits
122427ff
Commit
122427ff
authored
May 17, 2022
by
Tilman Aleman
Browse files
Merge branch 'main' of git.rwth-aachen.de:tilmandiego.aleman/reproduction
parents
80fbad5b
da964639
Changes
2
Hide whitespace changes
Inline
Side-by-side
helperfunctions.py
View file @
122427ff
...
...
@@ -3,21 +3,27 @@ from xfem.lsetcurv import *
import
matplotlib.pyplot
as
plt
from
netgen.csg
import
CSGeometry
,
OrthoBrick
,
Pnt
import
math
from
functools
import
reduce
def
refineMesh
(
mesh
,
phi
,
order
,
hs
=
10
**
6
):
if
order
>
2
:
hs
=
10
**
8
with
TaskManager
():
mesh
.
SetDeformation
(
None
)
lsetp1
=
GridFunction
(
H1
(
mesh
,
order
=
order
))
print
(
"lsetp1 interpolation.."
)
InterpolateToP1
(
phi
,
lsetp1
)
print
(
"lsetp1 refinement.."
)
RefineAtLevelSet
(
lsetp1
)
mesh
.
Refine
()
print
(
"lsetmeshadap..."
)
lsetmeshadap
=
LevelSetMeshAdaptation
(
mesh
,
order
=
order
,
threshold
=
1000
,
discontinuous_qn
=
True
,
heapsize
=
hs
)
print
(
"calc deform..."
)
deformation
=
lsetmeshadap
.
CalcDeformation
(
phi
)
mesh
.
SetDeformation
(
deformation
)
print
(
"set deform.."
)
if
order
>
1
:
mesh
.
SetDeformation
(
deformation
)
lset_approx
=
lsetmeshadap
.
lset_p1
return
lset_approx
,
deformation
...
...
@@ -116,12 +122,117 @@ def plot_errors(l2errs, l2errskf, geom, saveplot, order=2, refs=2, c=1):
plt
.
savefig
(
"convergence_plot_"
+
geom
+
"_c"
+
str
(
c
)
+
"_order"
+
str
(
order
)
+
"_refs"
+
str
(
refs
)
+
".png"
)
else
:
plt
.
show
()
def
computeEV
(
phi
,
mesh
,
lset_approx
,
deformation
,
order
=
2
):
vlams
=
[]
def
computeEV
(
phi
,
mesh
,
lset_approx
,
deformation
,
order
=
2
):
with
TaskManager
():
print
(
"starting ev comp"
)
vlams
=
[]
Vh
=
VectorH1
(
mesh
,
order
=
order
,
dirichlet
=
[],
low_order_space
=
True
)
ci
=
CutInfo
(
mesh
,
lset_approx
)
ba_IF
=
ci
.
GetElementsOfType
(
IF
)
VhG
=
Restrict
(
Vh
,
ba_IF
)
blocks
=
[]
freedofs
=
VhG
.
FreeDofs
()
class
SymmetricGS
(
BaseMatrix
):
def
__init__
(
self
,
smoother
):
super
(
SymmetricGS
,
self
).
__init__
()
self
.
smoother
=
smoother
def
Mult
(
self
,
x
,
y
):
y
[:]
=
0.0
self
.
smoother
.
Smooth
(
y
,
x
)
self
.
smoother
.
SmoothBack
(
y
,
x
)
def
Height
(
self
):
return
self
.
smoother
.
height
def
Width
(
self
):
return
self
.
smoother
.
height
vertexdofs
=
BitArray
(
Vh
.
ndof
)
vertexdofs
[:]
=
False
print
(
"doing stuff for precond"
)
"""
for v in mesh.vertices:
# vdofs = set(d for d in VhG.GetDofNrs(el) if freedofs[d])
for d in Vh.GetDofNrs(v):
vertexdofs[d] = True
# blocks.append(vdofs)
blocks = [set(d for el in mesh[v].elements for d in VhG.GetDofNrs(el) if freedofs[d]) for v in mesh.vertices]
# blocks = [reduce(lambda x,y: x or y,[set(d for d in VhG.GetDofNrs(el) if freedofs[d])for el in mesh[v].elements]) for v in mesh.vertices]
def blockcr(FEspace):
for v in FEspace.mesh.vectices:
vdofs = set()
for el in FEspace.mesh[v].elements:
vdofs |= set(d for d in FEspace.GetDofNrs(el)
if freedofs[d])
# blocks = pool.map(partial(calcstuff, mesh=mesh), range(len(mesh.vertices)))
vertexdofs &= VhG.FreeDofs()
"""
print
(
"done with stuff for precond"
)
n
=
Normalize
(
grad
(
lset_approx
))
ntilde
=
GridFunction
(
VhG
)
ntilde
.
Set
(
Normalize
(
grad
(
phi
,
mesh
)),
definedonelements
=
ba_IF
)
# n=ntilde
h
=
specialcf
.
mesh_size
# epsilon=1
eta
=
1
/
(
h
)
**
2
rhou
=
1.0
/
h
rhob
=
h
# Measure on surface
ds
=
dCut
(
lset_approx
,
IF
,
definedonelements
=
ba_IF
,
deformation
=
deformation
)
# Measure on the bulk around the surface
dX
=
dx
(
definedonelements
=
ba_IF
,
deformation
=
deformation
)
u2
=
VhG
.
TrialFunction
()
v2
=
VhG
.
TestFunction
()
# Weingarten mappings
ninter
=
GridFunction
(
VhG
)
ninter
.
Set
(
n
,
definedonelements
=
ba_IF
)
weing
=
grad
(
ninter
)
Pmat
=
Pmats
(
n
)
#.Compile()
def
E
(
u
):
return
((
Pmat
*
Sym
(
grad
(
u
))
*
Pmat
)
-
(
u
*
n
)
*
weing
)
#.Compile()
a2
=
BilinearForm
(
VhG
,
symmetric
=
True
)
a2
+=
InnerProduct
(
E
(
u2
),
E
(
v2
))
*
ds
a2
+=
(
Pmat
*
u2
)
*
(
Pmat
*
v2
)
*
ds
a2
+=
rhou
*
((
grad
(
u2
)
*
n
)
*
(
grad
(
v2
)
*
n
))
*
dX
a2
+=
(
eta
*
((
u2
*
ntilde
)
*
(
v2
*
ntilde
)))
*
ds
pre1
=
Preconditioner
(
a2
,
"direct"
,
inverse
=
"pardiso"
)
a2
.
Assemble
()
coarsepre
=
a2
.
mat
.
Inverse
(
vertexdofs
,
inverse
=
"pardiso"
)
# pre = Projector(VhG.FreeDofs(), True)
# pre1 = SymmetricGS(a2.mat.CreateSmoother(VhG.FreeDofs()))
# pre1 = SymmetricGS(a2.mat.CreateBlockSmoother(blocks, parallel=True))
pre
=
pre1
#+coarsepre
# pre = a2.mat.Inverse(freedofs=VhG.FreeDofs(), inverse="pardiso")
b
=
BilinearForm
(
VhG
,
symmetric
=
True
)
b
+=
InnerProduct
(
Pmat
*
u2
,
Pmat
*
v2
)
*
ds
b
+=
rhob
*
(
grad
(
u2
)
*
n
)
*
(
grad
(
v2
)
*
n
)
*
dX
b
+=
((
u2
*
n
)
*
(
v2
*
n
))
*
ds
b
.
Assemble
()
return
solvers
.
PINVIT
(
a2
.
mat
,
b
.
mat
,
pre
,
3
,
GramSchmidt
=
True
,
maxit
=
15
)
def
getevproblem
(
mesh
,
phi
,
lset_approx
,
deformation
,
order
=
2
):
Vh
=
VectorH1
(
mesh
,
order
=
order
,
dirichlet
=
[])
ci
=
CutInfo
(
mesh
,
lset_approx
)
ba_IF
=
ci
.
GetElementsOfType
(
IF
)
VhG
=
Restrict
(
Vh
,
ba_IF
)
...
...
@@ -158,13 +269,18 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2):
a2
+=
(
Pmat
*
u2
)
*
(
Pmat
*
v2
)
*
ds
a2
+=
rhou
*
((
grad
(
u2
)
*
n
)
*
(
grad
(
v2
)
*
n
))
*
dX
a2
+=
(
eta
*
((
u2
*
ntilde
)
*
(
v2
*
ntilde
)))
*
ds
pre
=
Preconditioner
(
a2
,
"
direct
"
)
a2
.
Assemble
()
pre
=
Preconditioner
(
a2
,
"
local
"
)
b
=
BilinearForm
(
VhG
,
symmetric
=
True
)
b
+=
InnerProduct
(
Pmat
*
u2
,
Pmat
*
v2
)
*
ds
b
+=
rhob
*
(
grad
(
u2
)
*
n
)
*
(
grad
(
v2
)
*
n
)
*
dX
b
+=
((
u2
*
n
)
*
(
v2
*
n
))
*
ds
return
a2
,
b
,
pre
,
Vh
,
VhG
def
computeEV2
(
Vh
,
VhG
,
a2
,
b
,
pre
):
Vh
.
Update
()
VhG
.
Update
()
print
(
"updated.."
)
a2
.
Assemble
()
b
.
Assemble
()
return
solvers
.
PINVIT
(
a2
.
mat
,
b
.
mat
,
pre
,
3
,
GramSchmidt
=
True
)
return
solvers
.
PINVIT
(
a2
.
mat
,
b
.
mat
,
pre
,
3
,
GramSchmidt
=
True
,
maxit
=
20
)
\ No newline at end of file
killing.py
View file @
122427ff
...
...
@@ -263,18 +263,22 @@ if __name__ == "__main__":
for
k
in
range
(
3
):
vlams
.
append
([])
for
i
in
range
(
n_cut_ref
+
1
):
print
(
"refinement is: "
,
i
)
if
i
==
0
:
# a,b,pre, Vh, VhG = getevproblem(mesh, phi, lset_approx, deformation, order)
pass
if
i
!=
0
:
lset_approx
,
deformation
=
refineMesh
(
mesh
,
phi
,
order
)
with
TaskManager
():
vlam
,
ev
=
computeEV
(
phi
,
mesh
,
lset_approx
,
deformation
,
order
)
vlam
,
ev
=
computeEV
(
phi
,
mesh
,
lset_approx
,
deformation
,
order
)
# vlam, ev = computeEV2(Vh, VhG, a, b, pre)
for
k
in
range
(
3
):
vlams
[
k
].
append
(
vlam
[
k
])
if
i
<
2
and
i
!=
0
:
symmetries
=
0
for
k
in
range
(
3
):
# if math.log(abs(1 - vlams[k][1]), 10) < math.log(abs(vlams[k][0] - vlams[k][1]), 10):
if
math
.
log
(
abs
(
1
-
vlams
[
k
][
1
]),
2
)
<-
(
order
+
1
)
*
(
i
+
1
):
if
math
.
log
(
abs
(
1
-
vlams
[
k
][
1
]),
2
)
<-
(
order
+
1
)
*
(
i
):
print
(
"corresponding log: "
)
print
(
math
.
log
(
abs
(
1
-
vlams
[
k
][
1
]),
2
))
print
(
"condition: "
,
-
(
order
+
1
)
*
(
i
+
1
))
...
...
@@ -287,12 +291,13 @@ if __name__ == "__main__":
temp
-=
(
vlams
[
j
][
-
1
]
-
vlams
[
j
][
-
2
])
**
2
/
(
(
vlams
[
j
][
-
1
]
-
vlams
[
j
][
-
2
])
-
(
vlams
[
j
][
-
2
]
-
vlams
[
j
][
-
3
]))
# if math.log(abs(1 - vlams[j][-1]), 10) < math.log(abs(temp - vlams[j][-1]), 10):
if
math
.
log
(
abs
(
1
-
temp
),
2
)
<
-
(
order
+
1
)
*
(
i
+
2
):
print
(
"corresponding log: "
)
print
(
math
.
log
(
abs
(
1
-
temp
),
2
))
print
(
"condition: "
,
-
(
order
+
1
)
*
(
i
+
1
))
print
(
"corresponding log: "
)
print
(
math
.
log
(
abs
(
1
-
temp
),
2
))
print
(
"condition: "
,
-
(
2
*
order
+
1
)
*
(
i
))
if
False
or
math
.
log
(
abs
(
1
-
temp
),
2
)
<
-
(
2
*
order
+
1
)
*
(
i
):
symmetries
+=
1
if
i
!=
0
:
if
True
and
i
!=
0
:
print
(
"starting solve.."
)
print
(
"symmetries: "
,
symmetries
)
with
TaskManager
():
...
...
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