Skip to content
Permalink
Browse files

Merge branch 'periodic_dg_implicit' into 'master'

assemble DG with periodic boundaries

See merge request !214
  • Loading branch information
JSchoeberl committed Aug 9, 2017
2 parents 0703290 + 1be0a58 commit fcc9a5dd0260a347c5c4bde1e6ea8f04de1b1989
Showing with 98 additions and 3 deletions.
  1. +48 −3 comp/bilinearform.cpp
  2. +50 −0 py_tutorials/DG/periodic.py
@@ -375,11 +375,24 @@ namespace ngcomp
{
//add dofs of neighbour elements as well
Array<DofId> dnums_dg;
Array<int> elnums_per;

for (int i = 0; i < nf; i++)
{
nbelems.SetSize(0);
ma->GetFacetElements(i,elnums);

// timerDG1.Stop();
if(elnums.Size() < 2)
{
int facet2 = ma->GetPeriodicFacet(i);
if(facet2 > i)
{
ma->GetFacetElements (facet2, elnums_per);
elnums.Append(elnums_per[0]);
}
}

for (int k=0; k<elnums.Size(); k++)
nbelems.Append(elnums[k]);

@@ -1229,6 +1242,8 @@ namespace ngcomp
{
LocalHeap lh = clh.Split();

Array<int> elnums_per(2, lh);

Array<int> dnums, dnums1, dnums2, elnums, fnums, vnums1, vnums2;
for (int i : r)
{
@@ -1239,8 +1254,24 @@ namespace ngcomp

ma->GetFacetElements(i,elnums);
el1 = elnums[0];

int fac2 = i;
// timerDG1.Stop();
if(elnums.Size() < 2)
{
int facet2 = ma->GetPeriodicFacet(i);
if(facet2 > i)
{
ma->GetFacetElements (facet2, elnums_per);
elnums.Append(elnums_per[0]);
fac2 = facet2;
}
else
continue;
}

if(elnums.Size()<2) continue;

el2 = elnums[1];

ElementId ei1(VOL, el1);
@@ -1250,7 +1281,7 @@ namespace ngcomp
int facnr1 = fnums.Pos(i);

fnums = ma->GetElFacets(ei2);
int facnr2 = fnums.Pos(i);
int facnr2 = fnums.Pos(fac2);

{
lock_guard<mutex> guard(printmatasstatus2_mutex);
@@ -1450,7 +1481,7 @@ namespace ngcomp
{
LocalHeap lh = clh.Split();

Array<int> dnums, dnums1, dnums2, elnums, fnums1, fnums2, vnums1, vnums2;
Array<int> dnums, dnums1, dnums2, elnums, elnums_per, fnums1, fnums2, vnums1, vnums2;
for (int el1 : r)
{
ElementId ei1(VOL, el1);
@@ -1460,6 +1491,20 @@ namespace ngcomp
HeapReset hr(lh);

ma->GetFacetElements(fnums1[facnr1],elnums);

int facet2 = fnums1[facnr1];
// timerDG1.Stop();
if(elnums.Size() < 2)
{
facet2 = ma->GetPeriodicFacet(fnums1[facnr1]);
if(facet2 != fnums1[facnr1])
{
ma->GetFacetElements (facet2, elnums_per);
elnums.Append(elnums_per[0]);
}
}


if (elnums.Size()<2)
{
ma->GetFacetSurfaceElements (fnums1[facnr1], elnums);
@@ -1532,7 +1577,7 @@ namespace ngcomp
int el2 = elnums[0] + elnums[1] - el1;
ElementId ei2(VOL, el2);
fnums2 = ma->GetElFacets(ei2);
int facnr2 = fnums2.Pos(fnums1[facnr1]);
int facnr2 = fnums2.Pos(facet2);

{
lock_guard<mutex> guard(printmatasstatus2_mutex);
@@ -0,0 +1,50 @@
from netgen.geom2d import *

periodic = SplineGeometry()
pnts = [ (0,0), (1,0), (1,1), (0,1) ]
pnums = [periodic.AppendPoint(*p) for p in pnts]

periodic.Append ( ["line", pnums[0], pnums[1]],bc="outer")
# This should be our master edge so we need to save its number.
lright = periodic.Append ( ["line", pnums[1], pnums[2]], bc="periodic")
periodic.Append ( ["line", pnums[2], pnums[3]], bc="outer")
# Slave boundaries must be defined in the same direction as master ones,
# this is why the the point numbers of this spline are defined in the reverse direction,
# leftdomain and rightdomain must therefore be switched as well!
# We use the master number as the copy argument to create a slave edge.
periodic.Append ( ["line", pnums[0], pnums[3]], leftdomain=0, rightdomain=1, copy=lright, bc="periodic")

from ngsolve import *
# from netgen.geom2d import unit_square

mesh = Mesh(periodic.GenerateMesh(maxh=0.05))


order=4
fes = L2(mesh, order=order, flags = { "dgjumps" : True })
u = fes.TrialFunction()
v = fes.TestFunction()

jump_u = u-u.Other()
jump_v = v-v.Other()
n = specialcf.normal(2)
mean_dudn = 0.5*n * (grad(u)+grad(u.Other()))
mean_dvdn = 0.5*n * (grad(v)+grad(v.Other()))

alpha = 4
h = specialcf.mesh_size
a = BilinearForm(fes)
a += SymbolicBFI(grad(u)*grad(v))
a += SymbolicBFI(alpha*order**2/h*jump_u*jump_v, skeleton=True)
a += SymbolicBFI(-mean_dudn*jump_v -mean_dvdn*jump_u, skeleton=True)
a += SymbolicBFI(alpha*order**2/h*u*v, BND, skeleton=True, definedon=mesh.Boundaries("outer"))
a += SymbolicBFI(-n*grad(u)*v-n*grad(v)*u, BND, skeleton=True, definedon=mesh.Boundaries("outer"))
a.Assemble()

f = LinearForm(fes)
f += SymbolicLFI(1*v)
f.Assemble()

gfu = GridFunction(fes, name="uDG")
gfu.vec.data = a.mat.Inverse() * f.vec
Draw (gfu)

0 comments on commit fcc9a5d

Please sign in to comment.
You can’t perform that action at this time.