

In [1]:

    
from mpi4py import MPI; import common; import grid; import functions; import la; import gdt
gdt.init_mpi()

$$\begin{align}\Omega &= [0, 1]^2\\ \Gamma_D &= \partial\Omega\end{align}$$



In [2]:

    
g = grid.make_cube_grid__2d_simplex_aluconform(lower_left=[0, 0], upper_right=[1, 1], num_elements=[4, 4], num_refinements=2, overlap_size=[0, 0])
#g.visualize('grid')

$$\begin{align}\kappa(x; \mu) &:= 1 - (1 - \mu) \cos(\tfrac{1}{2} \pi x_0) \cos(\tfrac{1}{2} \pi x_1)\\ f(x) &:= \tfrac{1}{2} \pi^2 \cos(\tfrac{1}{2} \pi x_0) \cos(\tfrac{1}{2} \pi x_1)\\ g_D(x) &:= 0\end{align}$$



In [3]:

    
#bump = functions.make_expression_function_1x1(g, 'x', 'cos(0.5*pi*x[0])*cos(0.5*pi*x[1])', order=3, name='bump')
#one = functions.make_constant_function_1x1(g, 1.0, name='one')
#diffusion = [one - bump, bump]
#diffusion[0].visualize(g, 'diffusion_affine_part')
#diffusion[1].visualize(g, 'diffusion_component')
#f = functions.make_expression_function_1x1(g, 'x', '0.5*pi*pi*cos(0.5*pi*x[0])*cos(0.5*pi*x[1])', order=3, name='rhs')
#f.visualize(g, 'force')
#g_D = functions.make_constant_function_1x1(g, 0.0, name='dirichlet')
kappa = functions.make_constant_function_1x1(g, 1.0, name='diffusion')
identity = functions.make_constant_function_2x2(g, [[0, 1], [1, 0]], name='id')
f = functions.make_constant_function_1x1(g, 1.0, name='force')
g_D = functions.make_constant_function_1x1(g, 0.0, name='dirichlet')
g_N = functions.make_constant_function_1x1(g, 0.0, name='neumann')



In [4]:

    
space = gdt.make_cg_space__1x1__p1__fem(g)
#space.visualize("cg_space")

elliptic_op = gdt.make_elliptic_matrix_operator__istl_sparse(kappa, space)

system_assembler = gdt.make_system_assembler(space)
system_assembler.append(elliptic_op)
system_assembler.assemble()