5/5 - (1 vote)

Part 1a: Horizontal Change

For and f(x) = 6x₁, the relationship

²u(x) = f(x),x (1)

can be shown to be true by taking the Laplacian of u(x).

The Python script used to discretize and solve the relationship is shown Appendix 1. The Paraview output of is shown in Figure 1.

The error reported by the script output is shown below:

errorL2 = 0.004931859322266561 errormax = 3.33066907388e16

Part 1b: Vertical Change

This relationship can be shown for and f(x) = 6x₂using the same steps as were used for and f(x) = 6x₁. The error from the python script used was the same as well. The Paraview output is shown in Figure 2, and the Python script used is given in Appendix 1b.

Figure 1: Paraview output for problem 1a

Figure 2: Paraview output for problem 1b

Problem 2: Improving Accuracy

Part 2a: Improving Through Mesh Resolution

The error was reduced by increasing the mesh resolution in the python script. Mesh resolutions of 4, 8, 16, 32, and 64 produced different errors. These are shown in Figure 3.

Figure 3: Error vs Mesh Sizing for Problem 2a

The code used to assess mesh refinement is provided in Appendix 2a.

Part 2b: Improving Through the Shape Function

The error was reduced by increasing the order of or the shape functions. Orders of 1, 2, 3, 4, and 5 were used. A plot of the error produced over the shape function order is given in Figure 4

The code used to assess shape function order is provided in Appendix 2b.

Figure 4: Error vs Mesh Sizing for Problem 2a

Appendix 1a: Code for Problem 1

Python script for 1a of Project 1

Original Author : Vinamra Agrawal

Date : January 25 , 2019

Edited By: Omkar Mulekar

Date : January 30 , 2019

Laplace (u) = f in the unit square u = u D on the boundary

In this problem , u D = x3

f = 6x

#============================================================

# This statement needs to be added for most ( i f not a l l )

# scripts for Fenics from future import print function from fenics import import matplotlib . pyplot as plt

#============================================================

#============================================================ # This is to create a mesh .

# In the class we saw how our rectangular domain was s p l i t # into triangles . This function , by defualt , creates a box # of dimensions 11 , with 8 intervals on each edge .

# These intervals are connected through triangles . mesh = UnitSquareMesh (8 , 8)

#============================================================

#============================================================ # This is to choose what kind of function space we are dealing # with . Essentially these are shape functions .

# Don t worry , we w i l l do shape functions in class later .

V = FunctionSpace (mesh , P , 1)

#============================================================

# This statement defines the boundary condition u D

# Expression is a convenient way to declare complicated # expressions in Fenics . On the backend , a program reads this # and converts i t in a format , the code can work with .

# This is a feature to make your l i f e easier . # degree = 2 t e l l s Fenics , the degree of the polynomial .

u D = Expression ( x [0] x [0] x [0] , degree=3)

#============================================================

# This is a function called boundary , that takes in two inputs

# x location of a point

# on boundary location of the boundary

# Right now, this function is not doing much. It is simply

# returning the location of the boundary . In the future , we w i l l # mess with this function as well . We w i l l see that i t is # extremely useful .

def boundary(x , on boundary ):

return on boundary

#============================================================

# This is where we actually define the boundary condition

# Fenics only needs you to define Dirichlet ( or displacement ) # boundary condition . You define the Neumann ( or force ) # boundary condition in a different way .

# This function takes in three inputs

# V the function space ( shape functions )

# u D the expression of the d i r i c h l e t boundary

# boundary a function ( declared above ) that t e l l s you

# exactly where the boundary is bc = DirichletBC (V, u D , boundary)

#============================================================

# This is where you	define	the weak form of the	equation
# First step : t e l l	Fenics	to start with u as a	t r i a l function
#u = TrialFunction (V)		in the function	space V
# Second step : t e l l	Fenics	to choose v(x) as a	test function
#		in the function	space V.

v = TestFunction (V)

# Third step : Define f (x) f = Expression ( 6x [0] , degree=1) # Fourth step : Define the l e f t hand side of the weak form a = dot( grad (u) , grad (v))dx

# Fifth step : Define the right hand side of the weak form L = fvdx

#============================================================

#============================================================ # This chunk of code computes the solution .

# First , move away from t r i a l function u and make i t a function . # Basically i t t e l l s Fenics to make u as a combination of # shape functions . We l l do this later too . u = Function (V)

# Finally . . . l e t i t run . Solve LHS = RHS, for u given BC. solve (a == L, u, bc)

#============================================================

#============================================================ # Save solution to f i l e in VTK format .

# In the class , we visualized the f i l e in Paraview .

vtkfile = File ( poisson/ solution . pvd ) vtkfile << u

#============================================================

# Now we test i f our code is working . This is where we use the

# method of manufactured solutions . We know that

# u = 1 + x2 + 2y2

# s a t i s f i e s the PDE. That s why we chose this as the boundary # condition . u = u D is automatically s a t i s f i e d .

# This line computes the L2 error between the computed solution

# u and the the solution we expect u D error L2 = errornorm (u D , u, L2 )

# Next , we find the value of the solution ( expected and computed) # at each vertex of the mesh . Then we compare the two . Our goal here # is to find the maximum error we get .

vertex values u D = u D. compute vertex values (mesh) vertex values u = u. compute vertex values (mesh)

import numpy as np

error max = np.max(np. abs( vertex values u D vertex values u ))

# Print errors print( error L2 = , errorL2 ) print( error max = , errormax )

Appendix 1b: Code for Problem 1

Python script for 1b of Project 1