poisson-pq2.py

Commented version of this example This example on gitlab

# SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
# SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
import numpy as np
import scipy.sparse.linalg
 
import dune.common
import dune.geometry
import dune.grid
import dune.functions
 
from dune.grid import DataType
from dune.functions import defaultGlobalBasis, Power, Composite, Lagrange
 
from fufemforms import *
from utilities import *
 
 
 
# Scalar right hand side
rhs = lambda x: 10
 
# Dirichlet boundary values
dirichletValues = lambda x : np.sin(2*np.pi*x[0])
 
def isNear(a,b):
  return np.abs(a-b) <= 1e-10
 
# Indicator function of Dirichlet boundary
isDirichlet = lambda x : isNear(x[0], 0) or isNear(x[0], 1) or isNear(x[1], 0) or isNear(x[1], 1)
 
 
 
# Number of grid elements (in one direction)
gridSize = 32
 
# Create a grid of the unit square
grid = dune.grid.structuredGrid([0,0],[1,1],[gridSize,gridSize])
 
 
 
# Create a nodal Lagrange FE basis
basis = defaultGlobalBasis(grid, Lagrange(order=2))
 
print("Dimension of FE space is "+str(len(basis)))
 
 
 
# Define Dirichlet constraints
N = len(basis)
 
# Mark all Dirichlet DOFs
isConstrained = np.zeros( N, dtype=bool )
basis.interpolate(isConstrained,isDirichlet)
 
# Interpolate the boundary values
constraintDOFValues = np.zeros( N )
basis.interpolate(constraintDOFValues, dirichletValues)
 
 
 
# Trial and test functions for variational problem
u = trialFunction(basis)
v = testFunction(basis)
 
# Scalar rhs coefficient function
f = Coefficient(rhs, basis.gridView)
 
# Bilinear form and rhs functional
a = integrate(dot(grad(u),grad(v)))
l = integrate(f*v)
 
 
 
# Assemble into matrix and vector
assembler = GlobalAssembler(basis)
A = assembler.assembleOperator(a)
b = assembler.assembleFunctional(l)
 
# Build constraints into matrix
incorporateEssentialConstraints(A, b, isConstrained, constraintDOFValues)
 
 
 
# Solve linear system!
x = scipy.sparse.linalg.spsolve(A, b)
 
 
u = basis.asFunction(x)
 
vtk = grid.vtkWriter(2)
u.addToVTKWriter("sol", vtk, DataType.PointData)
vtk.write("poisson-pq2-py")