The Poisson problem using Nitsche method on multiple subdomains

The Poisson problem using Nitsche method on multiple subdomains#

We consider a 2D domain \(\Omega\), that is subdivided into a grid of small squares, using the meshgrid function. Each subdomain has the form \((x_{i}, x_{i+1}) \times (y_{j}, y_{j+1})\), where \(x_1, ..., x_{n_x}\) and \(y_1, ..., y_{n_y}\) are subdivisions in each axis.

Formal Model#

from sympde.expr     import BilinearForm, LinearForm, integral, Norm
from sympde.expr     import find, EssentialBC
from sympde.topology import ScalarFunctionSpace, Square, Domain, element_of
from sympde.calculus import grad, dot
from sympde.calculus import jump, avg, minus, plus, Dn
from sympde.core     import Constant
from sympde.topology import meshgrid

from sympy import pi, sin

from psydac.api.discretization import discretize

# ... create a domain as meshgrid of a square
import numpy as np

x1 = np.linspace(0., 1., 3)
x2 = np.linspace(0., 1., 4)

domain = meshgrid(x1, x2)
# ...

# one sided approximation of the normal flux on the interface
Dn_I = lambda u: 0.5*(plus(Dn(u)) + minus(Dn(u)))

# internal interafaces of the domain
I = domain.interfaces

kappa = Constant('kappa', is_real=True)

V = ScalarFunctionSpace('V', domain)

x,y = domain.coordinates

u,v = [element_of(V, name=i) for i in ['u', 'v']]

# bilinear form
a = BilinearForm((u,v),
                   integral(domain, dot(grad(u),grad(v)))            
                 + integral(I, kappa * jump(u)*jump(v) - Dn_I(u)*jump(v) - jump(u)*Dn_I(v))
                )

# linear form
f = 2*pi**2*sin(pi*x)*sin(pi*y)
l = LinearForm(v, integral(domain, f*v))

# Dirichlet boundary conditions
bc = [EssentialBC(u, 0, domain.boundary)]

# Variational problem
equation   = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc)

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Cell In[1], line 7
      5 from sympde.calculus import jump, avg, minus, plus, Dn
      6 from sympde.core     import Constant
----> 7 from sympde.topology import meshgrid
      9 from sympy import pi, sin
     11 from psydac.api.discretization import discretize

ImportError: cannot import name 'meshgrid' from 'sympde.topology' (/opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/sympde/topology/__init__.py)

Discretization#

degree = [2,2]
ncells = [8,8]

# Create computational domain from topological domain
domain_h = discretize(domain, ncells=ncells, comm=None)

# Create discrete spline space
Vh = discretize(V, domain_h, degree=degree)

# Discretize equation
equation_h = discretize(equation, domain_h, [Vh, Vh])

Solving the PDE#

equation_h.set_solver('gmres', info=False, tol=1e-8)

uh = equation_h.solve(kappa=1e3)

Computing the Error Norm#

Computing the \(L^2\) norm#

ue = sin(pi*x)*sin(pi*y)

u = element_of(V, name='u')

# create the formal Norm object
l2norm = Norm(u - ue, domain, kind='l2')

# discretize the norm
l2norm_h = discretize(l2norm, domain_h, Vh)

# assemble the norm
l2_error = l2norm_h.assemble(u=uh)

# print the result
print(l2_error)
Contents