# Expression degree is not creating polynomial **URL:** \ **Category:** variational formulation\ **Created:** [March 4, 2021, 6:06pm UTC](https://fenicsproject.discourse.group/t/expression-degree-is-not-creating-polynomial/5224 "2021-03-04T18:06:41Z")\ **Posts on this page:** 11\ **Page:** 1
**Author:** ![mary\_h](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mary_h/32/1943_2.png) [@mary\_h](https://fenicsproject.discourse.group/u/mary_h)\ **Post date:** [March 4, 2021, 6:06pm UTC](https://fenicsproject.discourse.group/t/expression-degree-is-not-creating-polynomial/5224/1 "2021-03-04T18:06:41Z")
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?](https://fenicsproject.discourse.group/t/how-to-evaluate-a-userexpression-on-each-gauss-quadrature-point-during-assembly/3346)), 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 ``` ---
**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\ **Post date:** [March 4, 2021, 6:17pm UTC](https://fenicsproject.discourse.group/t/expression-degree-is-not-creating-polynomial/5224/2 "2021-03-04T18:17:53Z")
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?”?. ---
**Author:** ![mary\_h](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mary_h/32/1943_2.png) [@mary\_h](https://fenicsproject.discourse.group/u/mary_h)\ **Post date:** [March 4, 2021, 6:33pm UTC](https://fenicsproject.discourse.group/t/expression-degree-is-not-creating-polynomial/5224/3 "2021-03-04T18:33:42Z")
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, ---
**Author:** ![hernan](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/hernan/32/67_2.png) [@hernan](https://fenicsproject.discourse.group/u/hernan)\ **Post date:** [March 4, 2021, 11:09pm UTC](https://fenicsproject.discourse.group/t/expression-degree-is-not-creating-polynomial/5224/4 "2021-03-04T23:09:48Z")
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: ```python 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). ---
**Author:** ![mary\_h](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mary_h/32/1943_2.png) [@mary\_h](https://fenicsproject.discourse.group/u/mary_h)\ **Post date:** [March 4, 2021, 11:51pm UTC](https://fenicsproject.discourse.group/t/expression-degree-is-not-creating-polynomial/5224/5 "2021-03-04T23:51:39Z")
@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! ---
**Author:** ![hernan](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/hernan/32/67_2.png) [@hernan](https://fenicsproject.discourse.group/u/hernan)\ **Post date:** [March 4, 2021, 11:56pm UTC](https://fenicsproject.discourse.group/t/expression-degree-is-not-creating-polynomial/5224/6 "2021-03-04T23:56:42Z")
@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](https://fenicsproject.discourse.group/t/when-using-cg-degree-1-the-result-corresponds-to-degree-2/2751) for instance). ---
**Author:** ![mary\_h](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mary_h/32/1943_2.png) [@mary\_h](https://fenicsproject.discourse.group/u/mary_h)\ **Post date:** [March 5, 2021, 12:05am UTC](https://fenicsproject.discourse.group/t/expression-degree-is-not-creating-polynomial/5224/7 "2021-03-05T00:05:10Z")
@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! ---
**Author:** ![hernan](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/hernan/32/67_2.png) [@hernan](https://fenicsproject.discourse.group/u/hernan)\ **Post date:** [March 5, 2021, 12:06am UTC](https://fenicsproject.discourse.group/t/expression-degree-is-not-creating-polynomial/5224/8 "2021-03-05T00:06:44Z")
@mary_h I just edited my previous answer (you need to use the `write_checkpoint` function) ---
**Author:** ![mary\_h](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mary_h/32/1943_2.png) [@mary\_h](https://fenicsproject.discourse.group/u/mary_h)\ **Post date:** [March 5, 2021, 1:06am UTC](https://fenicsproject.discourse.group/t/expression-degree-is-not-creating-polynomial/5224/9 "2021-03-05T01:06:22Z")
@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! ---
**Author:** ![hernan](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/hernan/32/67_2.png) [@hernan](https://fenicsproject.discourse.group/u/hernan)\ **Post date:** [March 8, 2021, 1:59am UTC](https://fenicsproject.discourse.group/t/expression-degree-is-not-creating-polynomial/5224/10 "2021-03-08T01:59:48Z")
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. ---
**Author:** ![mary\_h](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mary_h/32/1943_2.png) [@mary\_h](https://fenicsproject.discourse.group/u/mary_h)\ **Post date:** [March 8, 2021, 2:17am UTC](https://fenicsproject.discourse.group/t/expression-degree-is-not-creating-polynomial/5224/11 "2021-03-08T02:17:51Z")
@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.