0
$\begingroup$

So I would like to make a custom color function that changes color depending on a threshold value, but also the amount of time the value of the function i am plotting, is above said threshold. I realize this might sound unclear so let me illustrate my problem :

Here is a code that solves a time dependant PDE, and plots its solution on a 3D cylinder.

r1 = 0.001;
r2 = 1;
l = 1;

tf = 200;

reg3D = ImplicitRegion[
   r1^2 <= x^2 + y^2 <= r2^2 && 0 <= z <= l, {x, y, z}];

eq = D[T[t, r, \[Theta], z], 
    t] - (1/r^2*D[T[t, r, \[Theta], z], \[Theta], \[Theta]] + 
     1/r*D[T[t, r, \[Theta], z], r] + D[T[t, r, \[Theta], z], r, r] + 
     D[T[t, r, \[Theta], z], z, z]);

sol = NDSolveValue[{eq == 
    NeumannValue[0, r == r1] + NeumannValue[0, r == r2], 
   T[0, r, \[Theta] , z] == 0, 
   DirichletCondition[T[t, r, \[Theta], z] == 0.6, 
    0 <= \[Theta] <= 2 Pi && r1 <= r <= 0.2], 
   PeriodicBoundaryCondition[
    T[t, r, \[Theta], z], \[Theta] == 2 Pi + 0.01 && 0 < z < l, 
    TranslationTransform[{0, -2 Pi, 0}]], 
   PeriodicBoundaryCondition[
    T[t, r, \[Theta], z], \[Theta] == -0.01 && 0 < z < l, 
    TranslationTransform[{0, 2 Pi, 0}]]}, 
  T, {r, r1, r2}, {\[Theta], -0.01, 2 *Pi + 0.01}, {z, 0, l}, {t, 0, 
   tf}]

Now you might notice that the equation i used here is rather trivial but that's just a toy model used to illustrate what I need.

Alright, on to the real aim of my question : I want to plot the result over a cylinder using the following criteria :

1 - if the local value of the solution is above a threshold value and thus, for a certain amount of time t0, the colorfunction returns Red

2 - Otherwise, the colorfunction returns White

Using the following code I managed to get something similar to what I wanted albeit with a simpler condition number 1, where I have disregarded the duration condition.

threshold = 0.5;
cf = If[# <= threshold, White, Red] &;
Table[SliceDensityPlot3D[
  sol[t, Sqrt[x^2 + y^2], Mod[ArcTan[x, y], 2 Pi], 
   z], {x^2 + y^2 == 0.98*r2^2, x^2 + y^2 == 1.01*r1^2, z == 0, 
   z == l}, {x, y, z} \[Element] reg3D, 
  ColorFunctionScaling  -> False, ColorFunction -> cf, Boxed -> False,
   Axes -> False, PlotPoints -> {50, 50, 25}, 
  ViewPoint -> Above], {t, {1, 2, 3, 4, 5, 10}}]  

enter image description here

That's great but as I have just written, this does not take into account the duration during which the function rises above the threshold, which needs to be above a certain duration t0 to be displayed Red.

So I have been experimenting for a week and a half now, trying to create the custom colorfunction I need and here is the best I could do through trials and error :

t0 = 10; 
threshold = 0.5;

ccf[t_, x_, y_] := Module[{t1, r1, \[Theta]1},
  t1 = t;
  Sow[x];
  Sow[y];
  r1 = Sqrt[x^2 + y^2];
  SetAttributes[sol, 
  HoldAll]; (*doesn't work without this line...for some reason*)
  \[Theta]1 = Mod[ArcTan[x, y], 2 Pi];
  Sow[t];
  Sow[r1];
  Sow[\[Theta]1];
  If[
   TrueQ[sol[t1, r1, \[Theta]1, l] >= threshold && t1 >= t0],
   Red,
   White
   ]
  ]

However this result in the following plots :

enter image description here

Now I've actually tried to understand why it was not working by putting some Sow[] and Reap[] to see what is being fed to my color function and from my beginner mathematica understanding, it seems that my time and space dependant colorfunction is being fed irrelevant values, or so that's how i interpret it.

One last point, the color function which I wrote seem to work outside of the plot code lines a point on the right of the cylinder is initially white and then red

but still isn't exactly what i needed. It simply compares the value of local variable t1 and check if it's above t0 which is different from checking the duration the function rises above the threshold as this would imply getting a the exact instant it first crosses that threshold and I have absolutely no idea how to do that. Besides, If I can't even make a simplified version work then there's no point in looking for something more complicated.

Therefore, I would like to ask you, how can i properly write the colorfunction that that i need. I would appreciate any help and insight on this matter.

Have a great day

$\endgroup$

1 Answer 1

1
+100
$\begingroup$

You have to embed the threshold information in the value of the function that you want to plot. Use a value that is far higher than the max of your solution in the whole domain, and add it to the solution when the condition is satisfied. You can then use a color function that takes this “bump” into account. In code :

threshold = 0.5;
thresholdIndicatorValue = 1000.;
t0 = 3.5
cf = If[# <=  threshold + thresholdIndicatorValue, 
    Opacity[0.1, White], Red] &;
decoratedSol = {t, x, y, z} |-> 
   sol[t, x, y, z] + thresholdIndicatorValue* Boole [t > t0];
Table[SliceDensityPlot3D[
  decoratedSol[t, Sqrt[x^2 + y^2], Mod[ArcTan[x, y], 2 Pi], 
   z], {x^2 + y^2 == 0.98*r2^2, x^2 + y^2 == 1.01*r1^2, z == 0, 
   z == l}, {x, y, z} \[Element] reg3D, ColorFunctionScaling -> False,
   ColorFunction -> cf, Boxed -> False, Axes -> False, 
  PlotPoints -> {50, 50, 25}, 
  ViewPoint -> Above], {t, {1, 2, 3, 4, 5, 10}}]

Regarding the condition that you want to achieve, it is equivalent to the following color function and decorated solution :

cf = If[# <=  thresholdIndicatorValue, Opacity[0.1, White], 
   Red] &;
decoratedSol = {t, x, y, z} |-> 
  sol[t, x, y, z] + 
   thresholdIndicatorValue* 
    Boole [NIntegrate[Boole[sol[s, x, y, z] > threshold], {s, 0, t}] >
       t0]; 

Unfortunately, you don’t want to use the code above because it will take an eternity to evaluate. The integration must be preprocessed at the level of the interpolating function instead.

This can be done the following way, using the utility packages for solution of NDSolve given by Wolfram :

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]
Needs["NDSolve`FEM`"]
reprocessNDSolveValueBySpaceSlice[thresholdValue_, solution_] := 
  With[
   {timeCoordinates = First@InterpolatingFunctionCoordinates[solution],
    timeLength = 
     First@Dimensions@InterpolatingFunctionValuesOnGrid[solution],
    spaceLength = 
     Last@Dimensions@InterpolatingFunctionValuesOnGrid[solution]
     },
   With[{reprocessingFunction = 
      Compile [{{values, _Real, 2}, {timePoints, _Real, 1}},
       Module[
        {i, j, 
         result = ConstantArray[0., {timeLength, spaceLength}]},
        For[i = 1, i < timeLength, i++,
         For[
          j = 1, j <= spaceLength, j++,
          result[[i + 1, j]] = result[[i, j]] +
            (timePoints[[i + 1]] - timePoints[[i]])*
             Piecewise[{
               {1,
                
                thresholdValue - values[[i, j]] < 0 && 
                 thresholdValue - values[[1 + i, j]] <= 
                  0}, {(-thresholdValue + 
                   values[[i, j]])/(values[[i, j]] - 
                   values[[1 + i, j]]),
                
                thresholdValue - values[[i, j]] < 0 && 
                 thresholdValue - values[[1 + i, j]] > 0},
               {(-thresholdValue + 
                   values[[1 + i, j]])/(-values[[i, j]] + 
                   values[[1 + i, j]]), 
                thresholdValue - values[[i, j]] >= 0 && 
                 thresholdValue - values[[1 + i, j]] < 0}},
              0]]
         ];
        result
        ]
       ]},
    ElementMeshInterpolation[(* build a new NDSolveValue-
     like Interpolation Function,*)
     InterpolatingFunctionCoordinates[
      solution],(*using the same coordinates in time and space,*)
     ArrayReshape[
      reprocessingFunction[(* and the reprocessing function, 
       applied to the values on the Grid*)
       InterpolatingFunctionValuesOnGrid[solution],
       timeCoordinates],
      {timeLength, 1, spaceLength}]
     ]
    ]
   ];
threshold = 0.5;
thresholdIndicatorValue = 1000.;
t0 = 3.5;
sol2 = reprocessNDSolveValueBySpaceSlice[threshold, sol];
cf = If[# < t0 &, White, Red] &;
Table[SliceDensityPlot3D[
  sol2[t, Sqrt[x^2 + y^2], Mod[ArcTan[x, y], 2 Pi], 
   z], {x^2 + y^2 == 0.98*r2^2, x^2 + y^2 == 1.01*r1^2, z == 0, 
   z == l}, {x, y, z} \[Element] reg3D, ColorFunctionScaling -> False,
   ColorFunction -> cf, Boxed -> False, Axes -> False, 
  PlotPoints -> {50, 50, 25}, 
  ViewPoint -> Above], {t, {5, 6, 7, 7.33, 7.66, 8}}]

One needs compilation because we are processing a 3D mesh, so a lot of points (around 38000 in your example), integrated over 74 time steps : the inner side of the loop is run millions of time ! The core of the computation is the Piecewise function. Indeed it is the result of an integration of

CriterionFunction = Boole[#>thresholdValue] &;

From 0 to 1, so the code can be generated by

Hold[Evaluate@
   Integrate[
    CriterionFunction[((value2 - value1) t + value1)], {t, 0, 
     1}]] /. {value1 -> values[[i, j]], 
  value2 -> ( values[[i + 1, j]])}

You could imagine a different criterion function to integrate symbolically and perform the same operation and replace the Piecewise function by the code that you obtained (without the Hold of course)

For instance, a “progressive criterion” like this


Piecewise[{{1, x > highThreshold}, 
  {(highThreshold - x)/(highThreshold - lowThreshold), 
   lowThreshold < x && x <= highThreshold}}, 0]

Would generate, using the code generation above (with a FullSimplify step after the Integration)

Hold[Piecewise[{{(2*(-highThreshold + values[[i,j]]) + 
      (highThreshold - values[[1 + i,j]])^2/(highThreshold - lowThreshold))/
     (2*(values[[i,j]] - values[[1 + i,j]])), 
    highThreshold < values[[i,j]] && highThreshold > values[[1 + i,j]] && 
     lowThreshold <= values[[1 + i,j]]}, 
   {(highThreshold + lowThreshold - 2*values[[i,j]])/
     (2*(-values[[i,j]] + values[[1 + i,j]])), 
    highThreshold < values[[i,j]] && lowThreshold > values[[1 + i,j]]}, 
   {1 + ((highThreshold - values[[i,j]])*(highThreshold - 2*lowThreshold + 
        values[[i,j]]))/(2*(highThreshold - lowThreshold)*
       (values[[i,j]] - values[[1 + i,j]])), lowThreshold < values[[i,j]] && 
     highThreshold < values[[1 + i,j]] && highThreshold > values[[i,j]]}, 
   {(highThreshold + lowThreshold - 2*values[[1 + i,j]])/
     (2*values[[i,j]] - 2*values[[1 + i,j]]), 
    lowThreshold >= values[[i,j]] && highThreshold < values[[1 + i,j]]}, 
   {1, (highThreshold == values[[i,j]] && values[[i,j]] < 
       values[[1 + i,j]]) || (highThreshold < values[[i,j]] && 
      highThreshold <= values[[1 + i,j]])}, 
   {((2*highThreshold - lowThreshold - values[[1 + i,j]])*
      (lowThreshold - values[[1 + i,j]]))/(2*(highThreshold - lowThreshold)*
      (values[[i,j]] - values[[1 + i,j]])), lowThreshold >= values[[i,j]] && 
     lowThreshold < values[[1 + i,j]] && highThreshold >= 
      values[[1 + i,j]]}, 
   {-1/2*(-2*highThreshold + values[[i,j]] + values[[1 + i,j]])/
      (highThreshold - lowThreshold), highThreshold >= values[[1 + i,j]] && 
     lowThreshold <= values[[1 + i,j]] && (highThreshold == values[[i,j]] || 
      (highThreshold >= values[[i,j]] && lowThreshold < values[[i,j]]))}, 
   {((2*highThreshold - lowThreshold - values[[i,j]])*
      (lowThreshold - values[[i,j]]))/(2*(highThreshold - lowThreshold)*
      (-values[[i,j]] + values[[1 + i,j]])), 
    lowThreshold > values[[1 + i,j]] && (highThreshold == values[[i,j]] || 
      (highThreshold >= values[[i,j]] && lowThreshold < values[[i,j]]))}}, 
  0]]

This is where symbolic integration shines…

Best

$\endgroup$
4
  • $\begingroup$ Thank you very much for you answer. I appreciate the opportunity to learn about this trick where you add a large value to create a sudden bump in the function values. The meaning in the condition I am trying to implement is, say we reason with temperatures for instance : I want to know if an area of the circle (or top of cylinder actually) has been "burnt" and my condition pretty much can be summed up by "if the temperature rises above X for too long, i.e for at least t0, then this area is burnt (in my case displayed as red). Perhaps this way of seeing it may be simpler to understand ? $\endgroup$ Feb 13, 2023 at 21:00
  • $\begingroup$ Ok. Give me 24/36 hours to code my idea… $\endgroup$ Feb 14, 2023 at 9:11
  • $\begingroup$ Hi, thank you for updating the answer you initially gave. I do have one problem however. When i use your first updated piece of code, i still obtain the constantly orange plots. I basically generated my sol and used your piece of code to generate and plot sol2 = reprocessNDSolveValueBySpaceSlice[threshold, sol];, which results in orange circles. Am I doing something wrong there ? $\endgroup$ Feb 21, 2023 at 9:04
  • $\begingroup$ After looking at your code in detail, i believe you should remove the first & in the colorfunction line of code cf = If[# < t0 &, White, Red] &;. Doing so plots the colors correctly. Thanks for the answer. $\endgroup$ Feb 22, 2023 at 12:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.