# Piecewise function to eliminate singularity but singularity still causing problem
**URL:** https://community.freefem.org/t/piecewise-function-to-eliminate-singularity-but-singularity-still-causing-problem/2235
**Category:** General Discussion
**Created:** [January 21, 2023, 4:05am UTC](https://community.freefem.org/t/piecewise-function-to-eliminate-singularity-but-singularity-still-causing-problem/2235 "2023-01-21T04:05:22Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@jsmacd](https://community.freefem.org/u/jsmacd)
#### Post date: [January 21, 2023, 4:05am UTC](https://community.freefem.org/t/piecewise-function-to-eliminate-singularity-but-singularity-still-causing-problem/2235/1 "2023-01-21T04:05:22Z")
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.
```auto
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;
```
---
### Author:  [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [January 25, 2023, 7:15am UTC](https://community.freefem.org/t/piecewise-function-to-eliminate-singularity-but-singularity-still-causing-problem/2235/2 "2023-01-25T07:15:32Z")
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.
---
### Author:  [@jsmacd](https://community.freefem.org/u/jsmacd)
#### Post date: [January 25, 2023, 3:25pm UTC](https://community.freefem.org/t/piecewise-function-to-eliminate-singularity-but-singularity-still-causing-problem/2235/3 "2023-01-25T15:25:36Z")
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++.
---
### Author:  [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [January 25, 2023, 3:43pm UTC](https://community.freefem.org/t/piecewise-function-to-eliminate-singularity-but-singularity-still-causing-problem/2235/4 "2023-01-25T15:43:52Z")
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.
---
### Author:  [@jsmacd](https://community.freefem.org/u/jsmacd)
#### Post date: [January 25, 2023, 9:15pm UTC](https://community.freefem.org/t/piecewise-function-to-eliminate-singularity-but-singularity-still-causing-problem/2235/5 "2023-01-25T21:15:12Z")
Now I understand how it is working! thank you 🙂