#! /usr/bin/env python3 # from dolfin import * def dpg_laplace ( my_grid, my_degree ): #*****************************************************************************80 # ## dpg_laplace applies DPG methods to the Poisson problem on the unit square. # # Discussion: # # This is a FEniCS implementation of a DPG method for the Dirichlet problem # # -div k grad u = f on UnitSquare # u = g on boundary, # # and is part of notes for graduate lectures introducing DPG methods. # # To my surprise, it seems that the degree of the Q space and element # must be set to 1. Using a value of 2 results in nan values. # JVB, 13 November 2018. # # In # Qs = FunctionSpace ( my_mesh, "BDM", 1, restriction = "facet" ) # I can replace the obscure 'restriction = "facet"' argument by # Qs = FunctionSpace ( my_mesh, "BDMF", 1 ) # and then I can use BDMF for the Qe statement as well. # JVB, 29 November 2018 # # I can access higher order BDMF elements if I impose a boundary condition. # I have tried zero boundary condition, and [-du/dx, -du/dy] condition. # However, the results are uniformly terrible. # JVB, 09 January 2019. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 09 January 2019 # # Modifier: # # John Burkardt # # Author: # # Jay Gopalakrishnan # # Reference: # # Jay Gopalakrishnan, # Five lectures on DPG Methods, # Spring 2013, Portland State University, # arXiv:1306.0557v2 [math.NA] 28 Aug 2014. # # Parameters: # # Input, integer MY_GRID, specifies the rows and columns of elements # in the mesh. # # Input, integer MY_DEGREE, specifies approximation degrees. # Error indicators: MY_DEGREE + 2 # Primary variable: MY_DEGREE + 1 # Fluxes: MY_DEGREE # import matplotlib.pyplot as plt # # Report input. # print ( '' ) print ( ' %dx%d mesh on the unit square' % ( my_grid, my_grid ) ) print ( " Approximation degree %d:" % ( my_degree ) ) # # Mesh the unit square. # my_mesh = UnitSquareMesh ( my_grid, my_grid, 'right/left' ) # # Plot the mesh. # label = ( 'dpg_laplace mesh (%d)' % ( my_grid ) ) plot ( my_mesh, title = label ) filename = ( 'dpg_laplace_mesh_%d.png' % ( my_grid ) ) plt.savefig ( filename ) print ( ' Graphics saved as "%s"' % ( filename ) ) plt.close ( ) # # Define the diffusivity function kappa(x). # kappa = Constant ( 1.0 ) # # Define the exact solution and right hand side. # exact = 1 if ( exact == 1 ): print ( ' Exact case #1' ) u_exact_expr = Expression ( "x[0]*(1-x[0])*x[1]*(1-x[1])", degree = 10 ) f_expr = Expression ( "2*x[1]*(1-x[1])+2*x[0]*(1-x[0])", degree = 10 ) elif ( exact == 2 ): print ( ' Exact case #2' ) u_exact_expr = Expression ( "exp(pow(x[0],2)+pow(x[1],2))*sin(x[0]+x[1])", degree = 10 ) f_expr = Expression ( "-2*exp(pow(x[0],2)+pow(x[1],2))*(sin(x[0]+x[1])*(1+2*(pow(x[0],2)+pow(x[1],2))) + 2*cos(x[0]+x[1])*(x[0]+x[1]))", degree = 10 ) else: print ( ' Exact case #3' ) u_exact_expr = Expression ( "2*(1+x[1])/(pow(3+x[0],2)+pow(1+x[1],2))", degree = 10 ) f_expr = Expression ( "0", degree = 0 ) # # Set spaces: # Es = Error estimator, # Us = primal variable, # Qs = interfacial flux. # Es = FunctionSpace ( my_mesh, "DG", my_degree + 2 ) Us = FunctionSpace ( my_mesh, "CG", my_degree + 1 ) Qs = FunctionSpace ( my_mesh, "BDMF", my_degree ) # # Set elements, needed for the MixedElement command: # Ee = Error estimator element, # Ue = primal variable element, # Qe = interfacial flux element. # Ee = FiniteElement ( "DG", triangle, my_degree + 2 ) Ue = FiniteElement ( "CG", triangle, my_degree + 1 ) Qe = FiniteElement ( "BDMF", triangle, my_degree ) # # Latest FENICS removed MixedFunctionSpace() command. # EUQe = MixedElement ( Ee, Ue, Qe ) EUQs = FunctionSpace ( my_mesh, EUQe ) # # Extract the individual trial and test function factors from the space. # Here, "e" for example is a symbol for typical trial functions from Es. # ( e, u, q ) = TrialFunctions ( EUQs ) ( et, ut, qt ) = TestFunctions ( EUQs ) # # Compute the normal vectors. # n = FacetNormal ( my_mesh ) # # Set the components of the saddle point problem: # # Define the inner product for the error estimator space. # yip = dot ( grad ( e ), grad ( et ) ) * dx + e * et * dx # # b( (u,q), et ) = ( kappa * grad u, grad et ) - < q.n, et > # This is an equation for u and q. # b1 = dot ( kappa * grad ( u ), grad ( et ) ) * dx \ - dot ( q ( '+' ), n ( '+' ) ) * ( et ( '+' ) - et ( '-' ) ) * dS \ - dot ( q, n ) * et * ds # # b( (ut,qt), e ) = ( kappa * grad ut, grad e ) - < qt.n, e > # This is an equation for e. # b2 = dot ( grad ( e ), kappa * grad ( ut ) ) * dx \ - ( e ( '+' ) - e ( '-' ) ) * dot ( qt ( '+' ), n ( '+' ) ) * dS \ - e * dot ( qt, n ) * ds # # Form the saddle point problem: # # yip + b1 = F * et * dx # b2 = 0 # a = yip + b1 + b2 b = f_expr * et * dx # # Define function G such that G \cdot n = g. # This is for case exact=1. # class BoundarySource(UserExpression): def __init__(self, mesh, **kwargs): self.mesh = mesh super().__init__(**kwargs) def eval_cell(self, values, x, ufc_cell): cell = Cell(self.mesh, ufc_cell.index) n = cell.normal(ufc_cell.local_facet) values[0] = - (x[1]+2*x[0]*x[1]*x[1]-x[1]*x[1]-2*x[0]*x[1])*n[0] values[1] = - (x[0]+2*x[0]*x[0]*x[1]-2*x[0]*x[1]-x[0]*x[0])*n[1] def value_shape(self): return (2,) G = BoundarySource ( my_mesh, degree = 2 ) # # Boundary conditions on CG and BDMF: # bc1 = DirichletBC ( EUQs.sub(1), u_exact_expr, DomainBoundary() ) # bc2 = DirichletBC ( EUQs.sub(2), G, DomainBoundary() ) # bcs = [ bc1, bc2 ] bcs = [ bc1 ] # # Solve. # euq = Function ( EUQs ) solve ( a == b, euq, bcs = bcs ) e, u, q = euq.split ( deepcopy = True ) # # Compute errors. # fmsh = refine ( refine ( my_mesh ) ) er_l2 = errornorm ( u_exact_expr, u, norm_type = 'L2', degree_rise = 3, mesh = fmsh ) er_h1 = errornorm ( u_exact_expr, u, norm_type = 'H1', degree_rise = 3, mesh = fmsh ) er_in = sqrt ( assemble ( ( dot ( grad ( e ), grad ( e ) ) + e * e ) * dx ) ) print ( " L2-norm of (u - uh) = %15.12f" % er_l2 ) print ( " H1-norm of (u - uh) = %15.12f" % er_h1 ) print ( " Error estimator norm = %15.12f" % er_in ) # # Compute the H1 projection of U (a standard Galerkin solution) # uu = TrialFunction ( Us ) vv = TestFunction ( Us ) aa = ( dot ( grad ( uu ), grad ( vv ) ) + uu * vv ) * dx bb = ( f_expr + u_exact_expr ) * vv * dx bc = DirichletBC ( Us, u_exact_expr, DomainBoundary() ) pu = Function ( Us ) solve ( aa == bb, pu, bcs = bc ) erp = errornorm ( u_exact_expr, pu, norm_type = 'H1', degree_rise = 3, mesh = fmsh ) print ( " Error in H1 Projection = %15.12f" % ( erp ) ) if ( 1.0e-12 < abs ( erp ) ): print ( " Error:Projection ratio = %15.12f" % ( er_h1 / erp ) ) print ( " H1-norm of (uh - proj) = %15.12f" % \ errornorm ( u, pu, norm_type = 'H1', degree_rise = 0 ) ) # # Terminate. # return def dpg_laplace_test ( ): #*****************************************************************************80 # ## dpg_laplace_test tests dpg_laplace. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 09 January 2019 # # Author: # # John Burkardt # import dolfin import platform import time print ( time.ctime ( time.time() ) ) # # Report level = only warnings or higher. # level = 30 set_log_level ( level ) print ( '' ) print ( 'dpg_laplace_test:' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' FENICS version %s'% ( dolfin.__version__ ) ) print ( " Discontinuous Petrov Galerkin method for the Poisson problem" ) print ( " on a unit square with 0 boundary conditions." ) for my_grid in ( 4, 8 ): for my_degree in range ( 1, 3 ): dpg_laplace ( my_grid, my_degree ) print ( '' ) print ( 'dpg_laplace_test:' ) print ( ' Normal end of execution.' ) print ( '' ) print ( time.ctime ( time.time() ) ) return if ( __name__ == '__main__' ): dpg_laplace_test ( )