Piecewise function to eliminate singularity but singularity still causing problem

General Discussion
1

Dear community,

I am trying to eliminate a singularity within a cube by replacing the function with the singularity with a piecewise function that has no singularity. However, I’m getting some funky errors including an attempt to divide by 0. But I’m trying very hard not to divide by zero. Can someone see where I’ve gone wrong?

Here is the code that’s giving me the error.

load "msh3" // to use cube
load "MUMPS" //

real myRadius = 0.3;
func uBC = (1/(((x - 0.5)^2 + (y - 0.5)^2)^0.5))*(((x - 0.5)^2 + (y - 0.5)^2)^0.5 > myRadius)
+ x*(((x - 0.5)^2 + (y - 0.5)^2)^0.5 <= myRadius);

// // loop on grid refinement
int[int] nN = [5, 10]; //, 20, 40];

// // initialise meshes
mesh3 Th;
mesh T2d;

// // loop through grid refinement
// //{
  for(int i = 0; i < nN.n; i++){
    // // construct meshes
    Th = cube(nN(i), nN(i), nN(i));

    fespace V3d(Th, P1);
    V3d uBCh = uBC;

    //plot(ueh3, Th, wait = true);
    plot(uBCh, Th, wait = true);

    // // 1: face y=0, 2: face x=1, 3: face y=1, 4: face x= 0, 5: face z=0, 6: face z=1
    T2d = square(nN(i), nN(i));

    fespace V2d(T2d, P1); // cuts
    V2d ucut;

    cout << "nn = " << nN(i) << endl;

    // take some cuts
    real[int] zcutVec = [0.2, 0.4, 0.6, 0.9];
    for(int j = 0; j< zcutVec.n; j++){
      ucut = uBC(x,y,zcutVec(j));
      plot(ucut, nbiso = 30, cmm=" numerical solution cut at z = " + zcutVec(j), wait = true, fill = true, value = true);
    }
  }
// // }

cout << "END OF PROGRAM" << endl;
2

Because you are dividing by zero in func uBC at any point where ((x - 0.5)^2 + (y - 0.5)^2)^0.5 = 0. If you change your mesh to one that does not have a node on the circle (for example int[int] nN = [7, 11];), you won’t divide by zero.

3

Thank you so much for your response!

I’m confused though because I thought that I defined a piecewise function that is linear inside the circle centred at (x,y) = (0.5, 0.5) with radius myRadius. I thought that this was how one defined piecewise functions in freeFEM++.

4

The way you are defining the piecewise function would be fine except that you are dividing by zero at the point (x,y) = (0.5, 0.5). If you remove this pole from your function you won’t have this issue.

5

Now I understand how it is working! thank you :slight_smile:

1 Like