# Building a structure out of cubes one block at a time

I have simplified code in which the motion of a point charge is simulated N times (one particle at a time) in a static electric field. When the point charge comes to rest (using WhenEvent boundary conditions), a Cube is set up at that location. However, I have designed this code such that when the cube is created at the resting position, a BoundaryMeshRegion structure is set up, and becomes a part of the boundary condition for the next simulated point charge. Doing so would allow you to, for instance, build a tower of blocks, one block at a time (if the starting point of the particle in each simulation is r, as seen in this image:

A problem arises sometimes when running the simulation. For example, when trying to make a 5 cube tower an error like this thrown:

NDSolve::nbnum1: The function value {0.700116,0.705879,-0.895985}\[Element]BoundaryMeshRegion[BooleanRegion[#1||#2||#3||#4||#5&,{Cube[{-5,-5,-5},0.03],Cube[{0.484458,0.482154,-0.115},0.03],Cube[{0.487374,0.491032,-0.130001},0.03],Cube[{0.487829,0.490812,-0.145002},0.03],Cube[{0.493027,0.497869,-0.160003},0.03]}]] is not True or False when the arguments are {0.,0.700116,-0.0000380001,0.705879,-0.0000183699,-0.895985,-0.0000352221,-0.0000380001,-524.417,-0.0000183699,<<3>>}.


and the blocks stop stacking properly. Strangely, I have watched this code work for 15 cubes and sometimes more. But it always seems to fail at 5 cubes, beyond 15 What can be done to eliminate this error, such that any number of cubes could be stacked (within bounds of the starting point of the simulation, obviously)?

Here is the code:

Clear["Global*"];
Needs["NDSolveFEM"];

q = 1.617733*10^-18;
r = BoundaryMeshRegion[
RegionDifference[Cuboid[{-4, -4, -6.5}, {4, 4, 0}],
Cuboid[{-4, -4, -0.1}, {4, 4, 0}]]];
pol = -1;

V0 = 5000;

(*build electric field*)
{U, {Fx, Fy, Fz}} =
NDSolveValue[{Laplacian[V[x, y, z], {x, y, z}] == 0,
DirichletCondition[
V[x, y, z] ==
pol*V0/2, (-0.1 <= z <= 0) && (-1 <= y <= 1) && (-1 <= x <= 1)],
DirichletCondition[
V[x, y, z] ==
0, (-6.5 <= z <= -6.4) && (-4 <= y <= 4) && (-4 <= x <=
4)]}, {V, -Grad[V[x, y, z], {x, y, z}]}, {x, y, z} \[Element]
r];

(*Visualization of vector field*)
v = Show[
VectorPlot3D[{Fx, Fy, Fz}, {x, -4, 4}, {y, -4, 4}, {z, -6.5, 0},
PlotTheme -> "Detailed", ColorFunction -> "Rainbow",
PerformanceGoal -> "Quality", VectorScale -> 0.05,
VectorPoints -> 7, PlotLegends -> Automatic],
RegionPlot3D[r, PlotStyle -> Opacity[0.1]]]

mass = 6.52*10^-11;(*particle mass*)
SeedRandom[1234];
n = 2; (*number of simulations/cube number*)
X0 = Take[RandomReal[{0.7, .71}, 2 n], n]; Y0 =
Take[RandomReal[{0.7, .71}, 2 n], -n]; Z0 =
RandomReal[{-.9, -0.89}, n];
rrv = RandomReal[{-0.0001, 0.0001}, 3 n];
u0 = Take[rrv, n]; v0 = Take[rrv, {n + 1, 2 n}];
w0 = Take[rrv, -n];

cubeTable = {Cube[{-5, -5, -5}, 0.03]};(*seed for cubeTable append*)
output = {};

(*Setup equations of motion and boundary conditions*)
SymplecticLeapfrog = {"SymplecticPartitionedRungeKutta",
"DifferenceOrder" -> 2, "PositionVariables" :> qvars}; time = {t, 0,
2}; tm = 10*10^3; partRad = 0.02; qvars = {X[t], Y[t], Z[t]};
eqn[i_] := {X''[t] == q/mass tm^2 Fx /. {x -> X[t], y -> Y[t],
z -> Z[t]},
Y''[t] == q/mass tm^2 Fy /. {x -> X[t], y -> Y[t], z -> Z[t]},
Z''[t] == (q/mass tm^2 Fz + 0.000002) /. {x -> X[t], y -> Y[t],
z -> Z[t]}, X[0] == X0[[i]], Y[0] == Y0[[i]], Z[0] == Z0[[i]],
X'[0] == u0[[i]], Y'[0] == v0[[i]], Z'[0] == w0[[i]],
WhenEvent[Z[t] >= -0.115, end = t; "StopIntegration"],
WhenEvent[{X[t], Y[t], Z[t]} \[Element] cubicMesh // Evaluate,
end = t; "StopIntegration"]};

h = 10^-4;

(*simulate the motion one cube at a time and build tower*)
Do[

cubicMesh = BoundaryMeshRegion[RegionUnion[cubeTable]];

soln[i] =
NDSolve[eqn[i], qvars, time, StartingStepSize -> h,
Method -> SymplecticLeapfrog];
T[i] = end;

AppendTo[output, qvars /. soln[i][[1]] /. t -> T[i]];
cubesMade = Cube[#, 0.03] & @output[[i]];

, {i, 1, n}]

(*show output*)
Show[Region[RegionUnion[cubeTable]], Graphics3D[{EdgeForm[None], r}],
Background -> LightGray, Boxed -> False,
PlotRange -> {{0, 1}, {0, 1}, {-1, 0}}]


• I'm not sure of the exact issue but it looks like it happens at BoundaryMeshRegion[RegionUnion[cubeTable] for certain values of cubeTable. It doesn't happen at n=5 for SeedRandom[1235], this makes me think there is something about the specific values in cubeTable when n=5 with SeedRandom[1234] (rather then the number of values in cubeTable). Perhaps it has issues if the cubes get too close to one another. Commented Aug 15, 2022 at 22:58
• @RudyPotter ya it is strange.
– Zach
Commented Aug 16, 2022 at 1:42
• Use ToElementMesh` and OpenCascadeLink for this, not BoundaryMeshRegion. Commented Aug 16, 2022 at 4:30
• you could try with cubicMesh = RegionUnion[BoundaryDiscretizeGraphics /@ cubeTable] Commented Aug 16, 2022 at 4:42