# Periodic boundary condition conflicts with DirichletCondition?

I want to solve a simple diffuse equation in a cylinder.

$$\nabla\cdot(\hat{c}\nabla\rho_P）=0$$

$$\rho_P(x,y,\theta)==0,\sqrt{x^2+y^2}=R$$

$$\rho_P(0,0,\theta)==1,\sqrt{x^2+y^2}=0$$

$$\rho_P(x,y,2\pi)=\rho_P(x,y,0)$$

$$\hat{c}=\begin{pmatrix} D_t& 0& 0 \\ 0 & D_t &0\\0 &0 &D_r \end{pmatrix} \quad$$

$$D_t: transportation\ coefficient$$

$$D_r: rotation\ coefficient$$

Needs["NDSolveFEM"]

{R, Difft, Diffr} = {29.5, 1, 0.01};
\[CapitalOmega] = Cylinder[{{0, 0, 0}, {0, 0, 2 \[Pi]}}, R];

c = {{Difft, 0, 0}, {0, Difft, 0}, {0, 0, Diffr}};

With[{rhop = rhop[x, y, \[Theta]]},
op1 = Div[c . Grad[rhop, {x, y, \[Theta]}], {x, y, \[Theta]}]];

Subscript[\[CapitalGamma], D1] =
DirichletCondition[rhop[x, y, \[Theta]] == 1,
x^2 + y^2 == 0 && \[Theta] != 2 \[Pi]];
Subscript[\[CapitalGamma], D2] =
DirichletCondition[rhop[x, y, \[Theta]] == 0,
x^2 + y^2 > 0 && \[Theta] != 0 && \[Theta] != 2 \[Pi]];
Subscript[\[CapitalGamma], p] =
PeriodicBoundaryCondition[rhop[x, y, \[Theta]], \[Theta] == 2 \[Pi],
Function[\[Theta], \[Theta] - 2 \[Pi]]];

cylinderMesh =
ToElementMesh[\[CapitalOmega], "MaxBoundaryCellMeasure" -> 0.04];

npfun = NDSolveValue[{op1 == 0, Subscript[\[CapitalGamma], p],
Subscript[\[CapitalGamma], D1], Subscript[\[CapitalGamma], D2]},
rhop, {x, y, \[Theta]} \[Element] cylinderMesh];

ContourPlot[npfun[x, y, 2 \[Pi]], {x, y} \[Element] Disk[{0, 0}, R],
PlotRange -> All]


Obviously, it is a symmetrical function, but when I solve it using NDSolveValue, the result is not symmetric.

I think I can solve it by adding

additionalPoints = {{0, 0, 0}, {0, 0, 2 \[Pi]}}
cylinderMesh =


But the error: 无法产生网络：enable to generate mesh

• Your second DirichletCondition should be Subscript[\[CapitalGamma], D2] = DirichletCondition[rhop[x, y, \[Theta]] == 0, x^2 + y^2 ==R^2 && \[Theta] != 0 && \[Theta] != 2 \[Pi]] Mar 16, 2023 at 9:38
• Region Cylinder suggest a 3D problem, but I'm missing an axial z coordinate!? Mar 16, 2023 at 11:05
• but,if i write as Subscript[\[CapitalGamma], D1] = DirichletCondition[rhop[x, y, \[Theta]] == 1, x^2 + y^2 == 0 && \[Theta] != 0 && \[Theta] != 2 \[Pi]];, Subscript[\[CapitalGamma], D2] = DirichletCondition[rhop[x, y, \[Theta]] == 0, x^2 + y^2 == R^2 && \[Theta] != 0 && \[Theta] != 2 \[Pi]]; there is a error: NDSolveValue::bcnop: No position was found on the boundary where x^2+y^2==0&&\[Theta]!=0&\[Theta]!=2 \[Pi] is true, so DirichletCondition[rhop==1,x^2+y^2==0&&\[Theta]!=0&&[Theta]!=2 \[Pi]]. will be ignored
– 江蛮子
Mar 17, 2023 at 0:58
• I think I should add a line(0,0,theta) or points(0,0,2pi)(0,0,0) to the mesh, but i can't find the option.
– 江蛮子
Mar 17, 2023 at 1:04
• The problem is not well-posed. Subscript[\[CapitalGamma], D2] already determines the solution, you cannot force b.c. at x^2+y^2==0, because it's the removable singularity in cylindrical coordinates. Related: mathematica.stackexchange.com/questions/245195/… Mar 17, 2023 at 12:44

I modified your code a little bit. You need two periodic boundaryconditions with correct TranslationTransform in all coordinates. Additionally I force NDSolve to use FiniteElement- Method (no need to load Needs["NDSolveFEM"])

I changed the Region \[CapitalOmega] to avoid singularity x^2+y^2->0. Play around with the new small parameter dR.

It is often recommended to rationalize the parameters.

{R, dR, Difft, Diffr} = {29.5, 29.5/100, 1, 0.01}// Rationalize[#, 0]&;
\[CapitalOmega] =RegionDifference[Cylinder[{{0, 0, 0}, {0, 0, 2 \[Pi] }},R], Cylinder[{{0, 0, 0}, {0, 0, 2 \[Pi] }}, dR]];


Try

c = {{Difft, 0, 0}, {0, Difft, 0}, {0, 0, Diffr}};

With[{rhop = rhop[x, y, \[Theta]]},op1 = Div[c . Grad[rhop, {x, y, \[Theta]}], {x, y, \[Theta]}]];

\[CapitalGamma]D1 =DirichletCondition[rhop[x, y, \[Theta]] == 1,x^2 + y^2 == dR^2 && 0 < \[Theta] < 2 \[Pi]];
\[CapitalGamma]D2 =DirichletCondition[rhop[x, y, \[Theta]] == 0,x^2 + y^2 == R^2 && 0 < \[Theta] < 2 \[Pi]];
\[CapitalGamma]p = {PeriodicBoundaryCondition[rhop[x, y, \[Theta]], \[Theta] >= 2 \[Pi] &&dR^2 < x^2 + y^2 < R^2, TranslationTransform[{0, 0, -2 Pi}]],PeriodicBoundaryCondition[rhop[x, y, \[Theta]], \[Theta] <= 0 && dR^2 <x^2 + y^2 < R^2, TranslationTransform[{0, 0, +2 Pi}]]};

npfun = NDSolveValue[{op1 ==0, \[CapitalGamma]p, \[CapitalGamma]D1, \[CapitalGamma]D2}, rhop, {x, y, \[Theta]} \[Element] \[CapitalOmega],
Method -> {"FiniteElement"}]

reg = ImplicitRegion[dR^2 <= x^2 + y^2 <= R^2 , {x, y}]
Plot3D[ npfun[x, y,  \[Pi]]  , {x, y} \[Element] reg,PlotRange -> All , Mesh -> False]


NDSolve now evaluates a solution which is plausible.

The periodic boundary conditions are fullfiled quite well:

Plot3D[Evaluate[npfun[x, y, 2 \[Pi]] - npfun[x, y, 0]], {x, y}\[Element] reg,PlotRange -> All]


Plot3D[Evaluate[Derivative[0, 0, 1][npfun][x, y, 2 \[Pi]] -Derivative[0, 0, 1][npfun][x, y, 0]], {x, y} \[Element] reg, PlotRange-> All]


I think you might get much better results if you formulate the underlying problem in "Cylindrical"-coordinates!