time series data

Featool multiphysics is really a good CFD software. While I can know the change of a variable with time from a plot in the GUI mode, may I extract the time series data in the .m code for further processing?
I tried using  c_data = { 'x', 'y', 'z','c','t'}; but it seems I can only extract the data at the end time point.

You should be able to use the evalexpr function in m-file scripts and on the Matlab command line to evaluate expressions in any point and time:

https://www.featool.com/doc/evalexpr_8m.html

Thank you for your reply. During simulations, can we enforce variables to be greater than or equal to zero (e.g. for variable such as concentration). Although we can do it in the post-processing steps, it should have been done during simulations with time and space. Thank you.

James wrote

Thank you for your reply. During simulations, can we enforce variables to be greater than or equal to zero (e.g. for variable such as concentration). Although we can do it in the post-processing steps, it should have been done during simulations with time and space. Thank you.

In the future, please start a new post for unrelated questions/issues.

As for your question, yes, if you try the equations and boundary conditions should accept the same syntax as for postprocessing expressions. However, the solver might not like it as your equation now is quite non-linear and potentially very unstable/oscillating.

I tried to calculate the variable 1 and variable_2 by calling my custom functions, where c is the concentration at all grid points of the test domain.
fea.expr = { 'Variable_1', 'myfun1(c)';'Variable_2', 'myfun2(c)'  };

While I can extract the time series concentration data at every grid point by using evalexpr (see below), can I get the time series data for Variable_1 and Variable_2 at all grid points? Now I have to re-run the custom functions after I have got the time-series data for concentration (c) at all grid points after the featool has run. Thank you!

c_data = { 'x', 'y'};
n_data = length(c_data);
p = fea.grid.p;
Positions = zeros(length(p),n_data);
for i=1:n_data
  Positions(:,i) = evalexpr( c_data{i}, p, fea );
end

You can get coordinates of the evaluation points by using a custom function for equation/boundary coefficinets of the form

  myfun( c, x, y )

Note that the coordinates will be in the quadrature points where the function is evaluated, and the function might be called multiple times for several blocks of cells in each time/evaluation step, as assembly and evaluation is performed vectorized over grouped blocks of grid cells.

If I have a simple test domain (say, Grid refined: 15 grid points, 16 triangulations) for a 2D diffusion/convection model.  
I used my custom function fea.expr = { 'Vxyt', 'myGridTest(c,x,y,t)'} to calculate the Vxyt from coordinates x, y and time t.
I found that in this custom function, there are c values at 16 (x,y) positions instead of 15 grid points.
i.e.
x1 y1 c1
....
x16 y16 c16
Does it mean these 16 positions represent 16 triangulations?
If each grid cell area is S and if I want to calculate the average c over the test domain in this custom function, I only need to use the following formula?

(c1*S + c2*S + ... + C16*S)/16*S
Thanks!

The evaluation functions are called during matrix assembly in the numerical quadrature points that are used to approximate integration. So it will depend on the quadrature rule you are using how many evaluation points there will be on each grid cell. In your example it looks like you might be using a first order rule with one evaluation point per cell. If all of your grid cells are the same size your formula should be correct in this specific case.