# Boundary conditions in a sphere using tetgreconstruction

**URL:** <https://community.freefem.org/t/boundary-conditions-in-a-sphere-using-tetgreconstruction/4043>\
**Category:** General Discussion\
**Created:** [August 21, 2025, 10:09am UTC](https://community.freefem.org/t/boundary-conditions-in-a-sphere-using-tetgreconstruction/4043 "2025-08-21T10:09:45Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![Joaquin](https://avatars.discourse-cdn.com/v4/letter/j/e9c0ed/32.png) [@Joaquin](https://community.freefem.org/u/Joaquin)\
**Post date:** [August 21, 2025, 10:09am UTC](https://community.freefem.org/t/boundary-conditions-in-a-sphere-using-tetgreconstruction/4043/1 "2025-08-21T10:09:45Z")

</div>

Hey, I have been trying to solve the Schrodinger equation for the hydrogen atom using Finite Elements in FreeFem. The problem is that in order to get good results I need to use a refine sphere because the Coulomb potential is stronger in the center of the sphere. When I am solving the variational equation I need to set some boundary conditions in the surface of the sphere but I do not know which boundary conditions I should set because when I build the refine sphere the surface does not get any labels and I am getting the next Warning : Warning: – Your set of boundary condition is incompatible with the mesh label. I do not know how to fix it if somebody could help me I appreciate it.

This is the code:

load “tetgen”

load “medit”

mesh Th=square(10,20,[x\*pi-pi/2,2\*y\*pi]); // a parametrization of a sphere  
func f1 =cos(x)\*cos(y);

func f2 =cos(x)\*sin(y);

func f3 = sin(x);

// partiel derivative of the parametrization DF

func f1x=sin(x)\*cos(y);  
func f1y=-cos(x)\*sin(y);

func f2x=-sin(x)\*sin(y);

func f2y=cos(x)\*cos(y);

func f3x=cos(x);

func f3y=0;

// M = DF^t DF

func m11=f1x^2+f2x^2+f3x^2;

func m21=f1x\*f1y+f2x\*f2y+f3x\*f3y;

func m22=f1y^2+f2y^2+f3y^2;

func perio=[[4,y],[2,y],[1,x],[3,x]];  
real hh=0.1;

real vv= 1/square(hh);

verbosity=2;

Th=adaptmesh(Th,m11\*vv,m21\*vv,m22\*vv,IsMetric=1,periodic=perio); Th=adaptmesh(Th,m11\*vv,m21\*vv,m22\*vv,IsMetric=1,periodic=perio);

plot(Th,wait=1);

verbosity=2;  
// construction of the surface of spheres

real Rmin = 1.;

func f1min = Rmin\*f1;

func f2min = Rmin\*f2;

func f3min = Rmin\*f3;  
meshS ThS=movemesh23(Th,transfo=[f1min,f2min,f3min]);

real[int] domain = [0.,0.,0.,145,0.01];  
mesh3 Th3sph=tetg(ThS,switch=“paAAQYY”,nbofregions=1,regionlist=domain);  
int[int] newlabel = [145,18];

real[int] domainrefine = [0.,0.,0.,145,0.0001];  
mesh3 Th3sphrefine=tetgreconstruction(Th3sph,switch=“raAQ”,region=newlabel,nbofregions=1,regionlist=domainrefine,sizeofvolume=0.0001);

int[int] newlabel2 = [145,53]; func fsize = 0.01/(( 1 + 5\*sqrt( (x-0.5)^2+(y-0.5)^2+(z-0.5)^2) )^3); mesh3 Th3sphrefine2=tetgreconstruction(Th3sph,switch=“raAQ”,region=newlabel2,sizeofvolume=fsize); medit(“sphere”,Th3sph,wait=1);

medit(“sphererefinedomain”,wait=1,Th3sphrefine);

medit(“sphererefinelocal”,wait=1,Th3sphrefine2);

int[int] L(Th3sphrefine2.nbe); // ahora L tiene tamaño = número de boundary elements for (int i=0; i\<Th3sphrefine2.nbe; i++) { L[i] = Th3sphrefine2.be(i).label; } cout \<\< “Boundary labels = " \<\< L \<\< endl; // --------------------------- // Constantes físicas // --------------------------- real me=9.10938356e-31; // masa electrón real planckh=6.62607004e-34; real hbar=1.0545718e-34; real kgm2Ds2ToeV=6.24150636309e18; real qe=1.6021766208e-19; // carga electrón real a0=5.2917721039e-11; // radio de Bohr real varepsilon0=8.854187817e-12; // permitividad // --------------------------- // Espacio de funciones // --------------------------- fespace Vh(Th3sphrefine,P2); Vh Psi, v; int nev=20; real tilE=-2e-18; Vh[int] EigenPsi(nev); real[int] EigenVal(nev); real R=20\*a0; // --------------------------- // Formulación variacional escalada // --------------------------- varf Schrodinger(Psi,v) = int3d(Th3sphrefine)( hbar^2/(2\*me\*R\*R)\*(dx(Psi)\*dx(v) + dy(Psi)\*dy(v) + dz(Psi)\*dz(v)) ) - int3d(Th3sphrefine)( qe^2/(4\*pi\*varepsilon0\*R\*sqrt(x^2+y^2+z^2))\*Psi\*v ) - int3d(Th3sphrefine)(tilE\*Psi\*v) + on(1,Psi=0); varf RHS(Psi,v) = int3d(Th3sphrefine)(Psi\*v); matrix A = Schrodinger(Vh,Vh,solver=“SPARSESOLVER”); matrix B = RHS(Vh,Vh,solver=CG,eps=1e-20); // --------------------------- // Resolución del problema // --------------------------- int num = EigenValue(A,B,sym=true,sigma=tilE, value=EigenVal,vector=EigenPsi, tol=1e-20,maxit=2000,ncv=0); // --------------------------- // Resultados // --------------------------- for(int i = 0; i \< nev; i++){ cout \<\< “Eigenvalue[”\<\<i\<\<”] = " \<\< EigenVal[i] \<\< " J , Energy level: " \<\< EigenVal[i]\*kgm2Ds2ToeV \<\< " eV" \<\< endl; }

---

<div class="post-metadata">

**Author:** ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)\
**Post date:** [August 25, 2025, 9:27am UTC](https://community.freefem.org/t/boundary-conditions-in-a-sphere-using-tetgreconstruction/4043/2 "2025-08-25T09:27:12Z")

</div>

Some corrections,

1. you make a error between label and region in 3d case.

I just put the label 1 on sphere (I do not understand why it is -1).

and the default solver say not enough memory so , I put MUMPS.

[aa-joaquin.edp](https://community.freefem.org/uploads/short-url/ttoYipt1V5lmkjiQQBlu9NfYcgc.edp) (3.8 KB)