# Balanced flux in FEA using NeumanValue

I'm using NeumannValue boundary conditions for a 3d FEA using NDSolveValue. In one area I have positive flux and in another area i have negative flux. In theory these should balance out (I set the flux inversely proportional to their relative areas) to a net flux of 0 but because of mesh and numerical inaccuracies they don't. Is there a way to constrain total flux = 0 and just set a constant flux for one of my areas?

edit: here's my boundary conditions:

Subscript[Γ, 1] =
NeumannValue[-1, (Abs[x] - 1)^2 + (Abs[y] - 1)^2 < (650/1000)^2 &&
z < -0.199  ];
Subscript[Γ, 2] =
NeumannValue[4, x^2 + y^2 + (z + 1/5)^2 < (650/1000/2)^2 ];


and my equations:

Dcof = 9000
ufun3d = NDSolveValue[
{D[u[t, x, y, z], t] - Dcof Laplacian[u[t, x, y, z], {x, y, z}] ==
Subscript[Γ, 1] + Subscript[Γ, 2],
u[0, x, y, z] == 0},
u, {t, 0, 10 }, {x, y, z} ∈ em];


and my element mesh:

a = ImplicitRegion[True, {{x, -1, 1}, {y, -1, 1}, {z, 0, 1}}];
b = Cylinder[{{0, 0, -1/5}, {0, 0, 0}}, (650/1000)/2];
c = Cylinder[{{1, 1, -1/5}, {1, 1, 0}}, 650/1000];
d = Cylinder[{{-1, 1, -1/5}, {-1, 1, 0}}, 650/1000];
e = Cylinder[{{1, -1, -1/5}, {1, -1, 0}}, 650/1000];
f = Cylinder[{{-1, -1, -1/5}, {-1, -1, 0}}, 650/1000];
r = RegionUnion[a,b,c,d,e,f];
boundingbox = ImplicitRegion[True, {{x, -1, 1}, {y, -1, 1}, {z, -1/5, 1}}];
r2 = RegionIntersection[r,boundingbox]
em = ToElementMesh[r2];


And this is what my mesh looks like from the bottom up. edit2: I figured i should add a plot of what i think is "wrong" too.
plotting the diagonal cross section i'd expect the values to be centered around 0 but they're all negative.

ContourPlot[ufun3d[5, xy, xy, z], {xy, -1 , 1 }, {z, -0.2, 1},
ClippingStyle -> Automatic, PlotLegends -> Automatic] • Can you help us to better help you by posting your code in a copy-and-paste-able manner? Jun 9, 2020 at 3:01
• I've updated with code. some variable names have been simplified. Jun 9, 2020 at 19:06
• @user1816847 I noticed that I was in error and deleted the comment, but note that your bottom cylinders have 16 times the area of the center cylinder. Shouldn't the factor 16 versus 4 if you want them to balance since flux is a per unit area quantity? Jun 9, 2020 at 23:55
• Indeed but each corner cylinder is only a 1/4 arc so its 16/4 = 4 times bigger than the center. Jun 10, 2020 at 0:26
• @user1816847 You need to use my user name so that I am informed of your reply. The code that you posted does not produce quarter arcs, rather it produces full cylinders. I posted an answer based on full cylinders. Jun 10, 2020 at 1:22

# Update (Steady-State Solution)

I think the fundamental issue is that you are over constraining your system. Whether you are solving the "heat equation" or not, your operator has the same form of the heat equation as shown below:

$$\rho {{\hat C}_p}\frac{{\partial T}}{{\partial t}} + \nabla \cdot {\mathbf{q}} = 0$$

If the flux, $$\mathbf{q}$$, needs to be perfectly conserved to conserve quanta, then it is equivalent to saying that the divergence of the flux is 0 or:

$$\nabla \cdot {\mathbf{q}} = 0$$

Therefore, the problem is a steady-state problem because there can be no accumulation in the domain:

$$\rho {{\hat C}_p}\frac{{\partial T}}{{\partial t}} + \nabla \cdot {\mathbf{q}} = \rho {{\hat C}_p}\frac{{\partial T}}{{\partial t}} + 0 = \rho {{\hat C}_p}\frac{{\partial T}}{{\partial t}} = 0$$

So, if you are seeing a response at all, then it is result of the numerical inaccuracies and not something physical.

If we substitute Fourier's Law for flux to put in terms of a temperature potential, we obtain:

$$\nabla \cdot {\mathbf{q}} = \nabla \cdot \left( { - {\mathbf{k}}\nabla T} \right) = \nabla \cdot \left( { - {\mathbf{k}}\nabla \left( {T + constant} \right)} \right)$$

The problem with this is that there is no unique solution because you can add an infinite number of constants to the temperature and still satisfy the equation. The way to obtain a unique solution is to add a Dirichlet or Robin condition on one of the boundaries and let the solver solve for the flux that balances the solution.

The following is a workflow that solves for the steady-state flux:

Needs["NDSolveFEM"]
Needs["OpenCascadeLink"]
a = ImplicitRegion[True, {{x, -1, 1}, {y, -1, 1}, {z, 0, 1}}];
b = Cylinder[{{0, 0, -1/5}, {0, 0, 0}}, (650/1000)/2];
c = Cylinder[{{1, 1, -1/5}, {1, 1, 0}}, 650/1000];
d = Cylinder[{{-1, 1, -1/5}, {-1, 1, 0}}, 650/1000];
e = Cylinder[{{1, -1, -1/5}, {1, -1, 0}}, 650/1000];
f = Cylinder[{{-1, -1, -1/5}, {-1, -1, 0}}, 650/1000];
shape0 = OpenCascadeShape[Cuboid[{-1, -1, 0}, {1, 1, 1}]];
shapeint = OpenCascadeShape[Cuboid[{-1, -1, -1}, {1, 1, 1}]];
union = OpenCascadeShapeUnion[shape0, shape1];
union = OpenCascadeShapeUnion[union, shape2];
union = OpenCascadeShapeUnion[union, shape3];
union = OpenCascadeShapeUnion[union, shape4];
union = OpenCascadeShapeUnion[union, shape5];
int = OpenCascadeShapeIntersection[union, shapeint];
groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp;
bmesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]
mesh = ToElementMesh[bmesh];
mesh["Wireframe"]
nv = NeumannValue[4, (x)^2 + (y)^2 < 1.01 (650/1000/2)^2 && z == -1/5];
dc = DirichletCondition[
u[x, y, z] == 0, (x)^2 + (y)^2 > 1.01 (650/1000/2)^2 && z == -1/5];
op = Inactive[
Div][{{-9000, 0, 0}, {0, -9000, 0}, {0, 0, -9000}}.Inactive[Grad][
u[x, y, z], {x, y, z}], {x, y, z}];
ufun3d = NDSolveValue[{op == nv, dc}, u, {x, y, z} \[Element] mesh];
ContourPlot[ufun3d[xy, xy, z], {xy, -Sqrt, Sqrt}, {z, -0.2, 1},
ClippingStyle -> Automatic, AspectRatio -> Automatic,
PlotLegends -> Automatic, PlotPoints -> {75, 50}]


The Mathematica (Top) result compares favorably to other FEM solver's such as Altair's AcuSolve (Bottom):

img = Uncompress[
"1:eJzt2+tP02cUB/\
CjYjQMnYuTYHQzLJItGI2OuWA0EpjG6eI07Vi8IFrgZ630Ai3VNjqeGQgCYyAKdlSBAuVS\
ZSgV5A5ekMWBEFEjYkBxBiUoTofxFvjamu2N/8GS8+KcnHOekzxvPm+\
Pb4ROtnMyERncaa1GoZR2TnS3Xq70vVEj6VWRwXq9whwxyTXwccUlV7hrPHyI3l50dKC5G\
ZWVKCpCdjYOHoTJhN27ERaGDRsQHIyAAPj5wccHnp4vp9Dwx9T3GXUtpvMrqeo7KtlMvyk\
peS/tSyTNYdpuI9nvtKqBvr5MX9ykOffJ8znRGw8a+YjuzqPuhdS6nGq+JcePdCyKfomj+\
AMUk0ERuRR6gtbU0rI2WnCdPh2gac8mTBifPv3p3Ll/+fvfCAz8Y/Xqerm8XKHIi41NF+\
LntDSD1SqVlm6qrl538eKKq1cX9ff7PnkyY2xsIkY/\
wOBs9HyOP5eiKQSnNiJPgUwtEvZjTwp2WbDVjvVOBJ3Dkk749mPmI0x+/\
WIqhrxxez6ufIlzQXCuR0E4sqKRZIY5CdFZCC/AxlMIacJX7Zh/G95DmPoCk8bg9RKz/\
sEnI/AbwqL7WNaH4B6suwZZJ7ZeRmQr1C0w1iO+\
CskVOORAjh0223hB3mjB8eFC673CnFtFRzuLslvtRxrtmc7iDEdJen5JmqU09dfS5MSyJH\
NZYowjQek4sO2ECK0Qm8+I7bVCahTRF4S+\
TZjaxU9dIuG6SOkRGX0ia0BYB4VtWJT8LcqfC+crUTsuml7HN4/ua35sbnqwt/\
GOsfGWoaE7tr5DV3dJU9cSXVunqnEqa8qls/\
aI6twdVZbwqkNhZ1K3OFPDKjMVFRblyXxNWbGhuNxU6Iy31SXktqRY29ItHVnZ3TmHe20Z\
A8VpD06mjJxOYk7MiTkxJ+\
bEnJgTc2JOzIk5MSfmxJyYE3NiTsyJOTEn5sScmBNzYk7MiTkxJ+\
bEnJgTc2JOzIk5MSfmxJyYE3NiTsyJOTEn5sScmBNzYk7MiTkxp/8dJ/\
kMIgrVGlRKrRS1VhsnKSV9oNzDNQwxx/17rOfuZEa1ZPB0Fd/\
o1Dq9PEYRKcndd3qyNSHvLX3436WfTDLo1MY4lU6rMrlm7625LwDd/+nVkmKPSqt89/\
KD3ii9BWHVFNA="];
dims = ImageDimensions[img];
colors2 =
RGBColor[#] & /@
ImageData[img][[IntegerPart@(dims[]/2), 1 ;; -1]];
DensityPlot[
ufun3d[X/Sqrt, X/Sqrt,
z], {X, -(Sqrt), (Sqrt)}, {z, -0.2, 1},
ColorFunction -> (Blend[colors2, #] &), PlotLegends -> Automatic,
PlotPoints -> {150, 100}, PlotRange -> All, AspectRatio -> Automatic,
Background -> Black, ImageSize -> Large] ## 3D Visualization Concepts

In the comments, @ABCDEMMM requested some 3D visualization of the solution. The example provided here, was actually quite complex as it appeared to have elements of clip-planes, iso-surfaces, and volume rendering. It is non-trivial to get all these elements tuned to produce a pleasing and informative visualization. In the process, I also could not get volume rendering (DensityPlot3D) and iso-surfaces (ContourPlot3D) to play nicely together. Here is an example workflow that combines clip-planes with volume rendering:

minmax = Chop@MinMax[ufun3d["ValuesOnGrid"]];
dpreg = DensityPlot3D[
ufun3d[x, y, z], {x, -1, 1}, {y, -1, 1}, {z, -0.2, 1},
PlotRange -> minmax, ColorFunction -> (Blend[colors2, #] &),
PlotLegends -> Automatic, OpacityFunction -> 0.05,
RegionFunction -> Function[{x, y, z, f}, -x + y > 0],
AspectRatio -> Automatic, Background -> Black, ImageSize -> Large]
dp = DensityPlot3D[
ufun3d[x, y, z], {x, -1, 1}, {y, -1, 1}, {z, -0.2, 1},
PlotRange -> minmax, ColorFunction -> (Blend[colors2, #] &),
PlotLegends -> Automatic, OpacityFunction -> 0.075,
AspectRatio -> Automatic, Background -> Black, ImageSize -> Large]
scp = SliceContourPlot3D[
ufun3d[x, y, z], {x == -0.9, y == 0.9, z == -0.15,
x - y == 0}, {x, -1, 1}, {y, -1, 1}, {z, -0.2, 1},
PlotRange -> minmax, Contours -> 30,
ColorFunction -> (Blend[colors2, #] &), PlotLegends -> Automatic,
RegionFunction -> Function[{x, y, z, f}, x - y <= 0.01],
AspectRatio -> Automatic, Background -> Black, ImageSize -> Large]
Show[dp, scp] Here is concept for 3D visualization using clip-planes and iso-surfaces:

cp100 = ContourPlot3D[
ufun3d[x, y, z], {x, -1, 1}, {y, -1, 1}, {z, -0.2, 1},
PlotRange -> minmax,
Contours -> (ufun3d[#/Sqrt, #/Sqrt, 0] & /@ {0.05, 0.32, 0.45,
0.65, 0.72, 0.78, 0.98}), MaxRecursion -> 0,
ColorFunctionScaling -> False,
ColorFunction -> (Directive[Opacity,
Blend[colors2, Rescale[#4, minmax]]] &), Mesh -> None,
PlotLegends -> Automatic, PlotPoints -> {100, 100, 50},
AspectRatio -> Automatic, Background -> Black, ImageSize -> Large]
cp50 = ContourPlot3D[
ufun3d[x, y, z], {x, -1, 1}, {y, -1, 1}, {z, -0.2, 1},
PlotRange -> minmax,
Contours -> (ufun3d[#/Sqrt, #/Sqrt, 0] & /@ {0.05, 0.32,
0.45, 0.65, 0.72, 0.78, 0.98}), MaxRecursion -> 0,
ColorFunctionScaling -> False,
ColorFunction -> (Directive[Opacity[0.5],
Blend[colors2, Rescale[#4, minmax]]] &), Mesh -> None,
PlotLegends -> Automatic, PlotPoints -> {100, 100, 50},
AspectRatio -> Automatic, Background -> Black, ImageSize -> Large];
cp25 = ContourPlot3D[
ufun3d[x, y, z], {x, -1, 1}, {y, -1, 1}, {z, -0.2, 1},
PlotRange -> minmax,
Contours -> (ufun3d[#/Sqrt, #/Sqrt, 0] & /@ {0.05, 0.32,
0.45, 0.65, 0.72, 0.78, 0.98}), MaxRecursion -> 0,
ColorFunctionScaling -> False,
ColorFunction -> (Directive[Opacity[0.25],
Blend[colors2, Rescale[#4, minmax]]] &), Mesh -> None,
PlotLegends -> Automatic, PlotPoints -> {100, 100, 50},
AspectRatio -> Automatic, Background -> Black, ImageSize -> Large];
scp25 = SliceContourPlot3D[
ufun3d[x, y, z], {x == -0.9, y == 0.9, z == -0.15, z == 0.90,
x - y == 0}, {x, -1, 1}, {y, -1, 1}, {z, -0.2, 1},
PlotRange -> minmax, Contours -> 30,
RegionFunction -> Function[{x, y, z, f}, x - y <= 0.1],
ColorFunction -> (Directive[Opacity[0.25], Blend[colors2, #]] &),
PlotLegends -> Automatic, PlotPoints -> {100, 100, 50},
AspectRatio -> Automatic, Background -> Black, ImageSize -> Large];
scp50 = SliceContourPlot3D[
ufun3d[x, y, z], {x == -0.9, y == 0.9, z == -0.15, z == 0.90,
x - y == 0}, {x, -1, 1}, {y, -1, 1}, {z, -0.2, 1},
PlotRange -> minmax, Contours -> 30,
RegionFunction -> Function[{x, y, z, f}, x - y <= 0.1],
ColorFunction -> (Directive[Opacity[0.5], Blend[colors2, #]] &),
PlotLegends -> Automatic, PlotPoints -> {100, 100, 50},
AspectRatio -> Automatic, Background -> Black, ImageSize -> Large];
scp100 = SliceContourPlot3D[
ufun3d[x, y, z], {x == -0.9, y == 0.9, z == -0.15, z == 0.90,
x - y == 0}, {x, -1, 1}, {y, -1, 1}, {z, -0.2, 1},
PlotRange -> minmax, Contours -> 30,
RegionFunction -> Function[{x, y, z, f}, x - y <= 0.1],
ColorFunction -> (Directive[Opacity, Blend[colors2, #]] &),
PlotLegends -> Automatic, PlotPoints -> {100, 100, 50},
AspectRatio -> Automatic, Background -> Black, ImageSize -> Large]
Show[scp50, cp25] It shows the 3D aspects of the solution and it is something to get you started. It will take time and practice to optimize the appearance of the plots.

# Update (Transient)

As alluded to in the comments, the $$t_{max} = 10$$ in the OP is about 18,000 times larger than it should be for a transient problem. One issue with running that long with a flux boundary condition is that the discretized areas of the boundary surfaces have an error associated with them that will accumulate with time. Therefore, one does not want to run more than necessary after the solution has reached a steady-state.

If we set the $$t_{max}=0.0001$$ and run the simulation with flux only boundary conditions, we can get a reasonable answer:

tmax = 0.0001;
nvin = NeumannValue[
4, (x)^2 + (y)^2 < 1.01 (650/1000/2)^2 && z == -1/5];
nvout = NeumannValue[-1, (x)^2 + (y)^2 > 1.01 (650/1000/2)^2 &&
z == -1/5];
ic = u[0, x, y, z] == 0;
op = Inactive[
Div][{{-9000, 0, 0}, {0, -9000, 0}, {0, 0, -9000}}.Inactive[Grad][
u[t, x, y, z], {x, y, z}], {x, y, z}] + D[u[t, x, y, z], t]
ufun3d = NDSolveValue[{op == nvin + nvout, ic},
u, {t, 0, tmax}, {x, y, z} ∈ mesh];
imgs = Rasterize[
DensityPlot[
ufun3d[#, X/Sqrt, X/Sqrt,
z], {X, -(Sqrt), (Sqrt)}, {z, -0.2, 1},
ColorFunction -> (Blend[colors2, #] &),
PlotLegends -> Automatic, PlotPoints -> {150, 100},
PlotRange -> All, AspectRatio -> Automatic, Background -> Black,
ImageSize -> Medium]] & /@ Subdivide[0, tmax, 30];
ListAnimate[imgs, ControlPlacement -> Top] As you can see, the density plot of the end point of the transient solution is essentially the same up to a constant as the previously calculated steady-state solution.

The code posted in the OP does not produce quarter arcs as suggested in the comments. On my machine, I obtain:

a = ImplicitRegion[True, {{x, -1, 1}, {y, -1, 1}, {z, 0, 1}}];
b = Cylinder[{{0, 0, -1/5}, {0, 0, 0}}, (650/1000)/2];
c = Cylinder[{{1, 1, -1/5}, {1, 1, 0}}, 650/1000];
d = Cylinder[{{-1, 1, -1/5}, {-1, 1, 0}}, 650/1000];
e = Cylinder[{{1, -1, -1/5}, {1, -1, 0}}, 650/1000];
f = Cylinder[{{-1, -1, -1/5}, {-1, -1, 0}}, 650/1000];
r = RegionUnion[a, b, c, d, e, f];
em = ToElementMesh[r];
em["Wireframe"] So, I am answering based on the full cylinders versus quarter arcs.

You will need a DirichletCondition or a Robin Condition somewhere to fully define temperature. Here is a case where applied a convective heat transfer condition to all but the bottom surfaces. There is a 16x change in area between the center port and the other ports, so I made the flux 16x more in the center. I also used the OpenCascadeLink to build the geometry since it seems to do a good job at snapping to features.

Needs["NDSolveFEM"]
a = ImplicitRegion[True, {{x, -1, 1}, {y, -1, 1}, {z, 0, 1}}];
b = Cylinder[{{0, 0, -1/5}, {0, 0, 0}}, (650/1000)/2];
c = Cylinder[{{1, 1, -1/5}, {1, 1, 0}}, 650/1000];
d = Cylinder[{{-1, 1, -1/5}, {-1, 1, 0}}, 650/1000];
e = Cylinder[{{1, -1, -1/5}, {1, -1, 0}}, 650/1000];
f = Cylinder[{{-1, -1, -1/5}, {-1, -1, 0}}, 650/1000];
shape0 = OpenCascadeShape[Cuboid[{-1, -1, 0}, {1, 1, 1}]];
union = OpenCascadeShapeUnion[shape0, shape1];
union = OpenCascadeShapeUnion[union, shape2];
union = OpenCascadeShapeUnion[union, shape3];
union = OpenCascadeShapeUnion[union, shape4];
union = OpenCascadeShapeUnion[union, shape5];
groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp;
bmesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]
mesh = ToElementMesh[bmesh];
mesh["Wireframe"]
nv1 = NeumannValue[-1/4, (x - 1)^2 + (y - 1)^2 < (650/1000)^2 &&
z < -0.199];
nv2 = NeumannValue[-1/4, (x + 1)^2 + (y - 1)^2 < (650/1000)^2 &&
z < -0.199];
nv3 = NeumannValue[-1/4, (x + 1)^2 + (y + 1)^2 < (650/1000)^2 &&
z < -0.199];
nv4 = NeumannValue[-1/4, (x - 1)^2 + (y + 1)^2 < (650/1000)^2 &&
z < -0.199];
nvc = NeumannValue[16,
x^2 + y^2 + (z + 1/5)^2 < (650/1000/2)^2 && z < -0.199];
nvconvective = NeumannValue[(0 - u[t, x, y, z]), z > -0.29];
ufun3d = NDSolveValue[{D[u[t, x, y, z], t] -
5 Laplacian[u[t, x, y, z], {x, y, z}] ==
nv1 + nv2 + nv3 + nv4 + nvc + nvconvective, u[0, x, y, z] == 0},
u, {t, 0, 10}, {x, y, z} \[Element] mesh];
ContourPlot[
ufun3d[5, xy, xy, z], {xy, -Sqrt, Sqrt}, {z, -0.2, 1},
ClippingStyle -> Automatic, PlotLegends -> Automatic,
PlotPoints -> 200] You could take advantage of symmetry and create 1/4 sized model. Here is a case where I applied a DirichletCondition to the top surface.

shaped = OpenCascadeShape[Cuboid[{0, 0, -1}, {2, 2, 2}]];
intersection = OpenCascadeShapeIntersection[union, shaped];
groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp;
bmesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]
mesh = ToElementMesh[bmesh];
mesh["Wireframe"]
nv1 = NeumannValue[-1/
4, (Abs[x] - 1)^2 + (Abs[y] - 1)^2 < (650/1000)^2 && z < -0.199];
nvc = NeumannValue[16/4,
x^2 + y^2 + (z + 1/5)^2 < (650/1000/2)^2 && z < -0.199];
dc = DirichletCondition[u[t, x, y, z] == 0, z == 1];
ufun3d = NDSolveValue[{D[u[t, x, y, z], t] -
5 Laplacian[u[t, x, y, z], {x, y, z}] == nv1 + nvc , dc,
u[0, x, y, z] == 0}, u, {t, 0, 10}, {x, y, z} ∈ mesh];
ContourPlot[ufun3d[5, xy, xy, z], {xy, 0, Sqrt}, {z, -0.2, 1},
ClippingStyle -> Automatic, PlotLegends -> Automatic] • I missed the region intersect when changing variablenames sorry. In my case it's not a heat but another quantity and the physical system makes sure that each quanta generated is also removed. The renaming boundary conditions also have 0 flux because of the physics. Jun 10, 2020 at 3:31
• also i didn't know about OpenCascadeShape so thanks for that! Jun 10, 2020 at 3:57
• @user1816847 I came across a potential issue that I am attempting to resolve. Normally, I get pretty good agreement between Mathematica and other solvers, such as COMSOL, but I have not done any 3D cases. I had to multiply the outer cylinder fluxes by 4x to get them to agree. This could be causing your balance issue. Also, do you have a value of Dcof? It is not defined in your code. Jun 11, 2020 at 1:17
• Dcof is 9000 (very fast). For me i have to set my fluxes to be -1 and about 4.7 to 4.8 which seems too far for numerical error too :/ Jun 11, 2020 at 2:57
• @user1816847 Please check out the updates. As discussed in the steady-state update, the flux is balanced by the solver using a flux condition and a Dirichlet condition. The solution is unique up to an arbitrary constant. As I showed, you can solve the transient problem too, but you should not run it much longer than necessary to evolve to steady-state. Otherwise, just the discretization error in the areas of the inflow and outflow will accumulate. Jun 13, 2020 at 0:46

Too long for a comment. An easy way to generate a high quality mesh is to replace the Implicitegion with Cubuid and make use of the OpenCascade boundary mesh generator:

Needs["NDSolveFEM"]
(*a=ImplicitRegion[True,{{x,-1,1},{y,-1,1},{z,0,1}}];*)

a = Cuboid[{-1, -1, 0}, {1, 1, 1}];
b = Cylinder[{{0, 0, -1/5}, {0, 0, 0}}, (650/1000)/2];
c = Cylinder[{{1, 1, -1/5}, {1, 1, 0}}, 650/1000];
d = Cylinder[{{-1, 1, -1/5}, {-1, 1, 0}}, 650/1000];
e = Cylinder[{{1, -1, -1/5}, {1, -1, 0}}, 650/1000];
f = Cylinder[{{-1, -1, -1/5}, {-1, -1, 0}}, 650/1000];
r = RegionUnion[a, b, c, d, e, f];
(*boundingbox=ImplicitRegion[True,{{x,-1,1},{y,-1,1},{z,-1/5,1}}];*)

boundingbox = Cuboid[{-1, -1, -1}, {1, 1, 1}];
r2 = RegionIntersection[r, boundingbox];

mesh = ToElementMesh[r2, "BoundaryMeshGenerator" -> {"OpenCascade"}];
groups = mesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp;
mesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]] We can use mesh of first order for 3D visualization and short time for visibility. We also change boundary conditions:

Needs["NDSolveFEM"]; a =
ImplicitRegion[True, {{x, -1, 1}, {y, -1, 1}, {z, 0, 1}}];
b = Cylinder[{{0, 0, -1/5}, {0, 0, 0}}, (650/1000)/2];
c = Cylinder[{{1, 1, -1/5}, {1, 1, 0}}, 650/1000];
d = Cylinder[{{-1, 1, -1/5}, {-1, 1, 0}}, 650/1000];
e = Cylinder[{{1, -1, -1/5}, {1, -1, 0}}, 650/1000];
f = Cylinder[{{-1, -1, -1/5}, {-1, -1, 0}}, 650/1000];
r = RegionUnion[a, b, c, d, e, f];
boundingbox =
ImplicitRegion[True, {{x, -1, 1}, {y, -1, 1}, {z, -1/5, 1}}];
r2 = RegionIntersection[r, boundingbox];
em = ToElementMesh[r2, "MeshOrder" -> 1, MaxCellMeasure -> 10^-4];
Subscript[\[CapitalGamma], 1] =
NeumannValue[-1, z == -1/5 && x^2 + y^2 > (650/1000/2)^2];
Subscript[\[CapitalGamma], 2] =
NeumannValue[4, z == -1/5 && x^2 + y^2 < (650/1000/2)^2]; Dcof = 9000;
ufun3d = NDSolveValue[{D[u[t, x, y, z], t] -
Dcof Laplacian[u[t, x, y, z], {x, y, z}] ==
Subscript[\[CapitalGamma], 1] + Subscript[\[CapitalGamma], 2],
u[0, x, y, z] == 0}, u, {t, 0, 10^-3}, {x, y, z} \[Element] em];

DensityPlot3D[
ufun3d[1/1000, x, y, z], {x, 0, 1}, {y, 0, 1}, {z, -1, 1},
ColorFunction -> "Rainbow", OpacityFunction -> None,
BoxRatios -> {1, 1, 1}, PlotPoints -> 50, Boxed -> False,
PlotLegends -> Automatic, Axes -> False] General view of 3D distribution from different points

DensityPlot3D[ufun3d[1/1000, x, y, z], {x, y, z} \[Element] em,
ColorFunction -> "Rainbow", OpacityFunction -> None,
BoxRatios -> Automatic, PlotPoints -> 50, Boxed -> False,
Axes -> False] 