Expression degree is not creating polynomial

variational formulation
1

Hello everyone,

I am trying to fit a cubic polynomial into my step function, and I have encountered a problem with expression degree. In my attached minimal code, I have an elliptical crack and for my crack vector, I am assigning the value of one when inside the crack, and zero for outside the crack.
My problem is when I use expression degree = 1 (d_initial= EllipticalCrackExpression(degree=1)), I get a step function with sharp endpoints. However, when I change the degree of expression to 3 in order to get a cubic polynomial, I am getting some oscillating crack vector with sharp points. Now, I am wondering how I can get a cubic curve that fits into my step function?
I have seen this post (How to evaluate a UserExpression on each gauss-quadrature point during assembly?), which I think has the same issue as mine. Here, as a solution, it is mentioned that we can use quadrature elements, but since I will use a variational problem, I think it might be buggy based on what has been mentioned in the response.

Thanks in advance for your great help and support

import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
parameters["form_compiler"]["quadrature_degree"] = 5

mesh = RectangleMesh(Point(-1,-1),Point(1,1), 100, 100)
V = VectorFunctionSpace(mesh, "CG", 1)

################################################################################
# Expression for elliptical crack geometry
################################################################################
class EllipticalCrackExpression(UserExpression):
    def eval(self, value, x):
        
        a=0.5
        b=0.1
        A = 1.2
        a_1 = 2/((A-1)**3)
        b_1 =-(3/2) *(1+A) * a_1 
        c_1= -3*a_1 -2*b_1
        d_1 = 1 - a_1 - b_1 - c_1

        if (x[0]/a)**2 + (x[1]/b)**2 <=1:
            value[0] = 0
            value[1] = 1 
        elif (x[0]/a)**2 + (x[1]/b)**2 >= A:
            value[0] = 0
            value[1] = 0 
        else:
            value[0] = 0 
            value[1] = (a_1 *((x[0] /a)**2 + (x[1]/b)**2)**3 \
                          + b_1*((x[0] /a)**2 + (x[1]/b)**2)**2 \
                              + c_1* ((x[0] /a)**2 + ( x[1]/b)**2) + d_1)
            
            
    def value_shape(self):
        return (2,)

d = Function(V)
d_initial= EllipticalCrackExpression(degree=1) #degree=3
d=project(d_initial,V)

d.rename("d", "crack_vector")
fileD = File("Crack_vector.pvd")
fileD << d
2

Apologies for being lazy and not reading your post in detail.

Is your question essentially “how can I avoid Gibbs phenomenon when interpolating data to finite element spaces?”?.

3

Many thanks for your response. Yes, that’s my question. Also, please note that I will use a variational problem in the rest of my code.

Thanks again,

4

The project function internally minimizes the L2 norm between a function in V and your UserExpression, i.e., solves the linear problem of finding w\in V such that

\int_{\Omega} u \cdot v ~ d\Omega = \int_{\Omega} f \cdot v ~ d\Omega, \qquad \forall v \in V,

which I believe is the cause of your issue.

To avoid that you can simply interpolate your expression into V using:

d=interpolate(d_initial,V)

which evaluates your expression on the degrees of freedom of V. However, I cannot understand why you want to project a UserExpression of higher degree into a CG1 space (the accuracy that you get using degree=3 is lost with the projection/interpolation).

1 Like
5

@hernan, thank you very much for your response. You are right, I should use a CG3 space. However, after I tried interpolation and space of CG3, although I am not getting the Gibbs phenomenon, I am still not getting the crack vector as a cubic curve, and its shape is very close to a degree =1 step function. Other than plotting the crack vector and visually checking that, is there any other way to make sure that the vector d is a cubic curve?

Thanks a lot for your help!

6

@mary_h I believe that when exporting data to vtk files, dolfin takes the nodal dofs of the output if the function belongs to a CG space.

To get a function for visualization with element-wise P3 polynomials you should try the write_checkpoint method of the XDMFFile class (see the docken answer in this post for instance).

7

@hernan thanks for your reply. I also tried with XDMF format, but still, I am not getting a cubic curve for the vector. I am viewing the file using Paraview, could that be a problem? If yes, can you please tell me how you are viewing your output file?

Thanks so much!

8

@mary_h I just edited my previous answer (you need to use the write_checkpoint function)

9

@hernan Thanks for your responses. Now, I am using this to save the file:

fileD = XDMFFile(“Crack_vector.xdmf”)
fileD.write_checkpoint(d, “d”, 2)
fileD.close()

But, I still don’t get the cubic curve. I even tried replacing VectorFunctionSpace by V = VectorFunctionSpace(mesh, "DG", 3) , but it is not working neither.

Is there anything else that I should try?
Thanks!

10

I don’t know why isn’t working. As far as I have read, write_checkpoint should work to visualize CG and DG functions of orders 1,2, and 3. Maybe the developers can give you some insights.

11

@hernan thanks for your response. I think my problem is solved. I had to try with a smaller domain to see the curve.
Thanks again for your great support and response.