Question about Compositional fields

ASPECT
#1

Hi all,

Recently i wrote a postprocess/visualization plugin, an important step to make this plugin work is to select things like material properties for different compositional fields.
Supposing there are serveral compositional fields named: stress_xx, stress_yy, stress_xy, upper_crust, lower_crust, inclusion, lithosphere.
Here is how i did in the function: evaluate_vector_field()
…
for(unsigned int q=0; q<n_quadrature_points; ++q)
{
if(this->introspection().name_for_compositional_index(q) == “upper_crust”)
{
statement…
}
else if(this->introspection().name_for_compositional_index(q) == “lower_crust”)
{
statement…
}
else if(this->introspection().name_for_compositional_index(q) == “lithosphere”)
{
statement…
}
}
However, an error occurred:
Index 7 is not in the half-open range [0,7)
How should I deal with this problem?Or is there a better way to realize such function?

1617438777833
1617438777833998×334 148 KB

Best,
Xie

#2

This is a semantic error. ‘q’ is an index that loops over quadrature points,
but the function name_for_compositional_index() requires an argument that
denotes an index that corresponds to a component of the solution vector. These
two don’t go together, and the error is a result of that.

Best
Wolfgang

#3

Thanks for replying. You’re right. I’d like to ask further. I found that I inevitably used nested loops, but this brought new troubles, as in the previous post:

for(unsigned int q=0; q<n_quadrature_points; ++q)
{
  for(unsigned int i=0; i<this->introspection().n_compositional_fields(); ++i)
  {
    if(this->introspection().name_for_compositional_index(q) == “upper_crust”)
    {
      statement...
    }
    else if(this->introspection().name_for_compositional_index(q) == “lower_crust”)
    {
      statement...
    }
    else if(this->introspection().name_for_compositional_index(q) == “lithosphere”)
    {
      statement...
    }
  }
  statement...
}

But it is obvious that the result of the calculation will be repeated many times, because of the number of cycles. I also tried other methods, but failed to get the results i wanted.
From my perspective, the key point of the issue seems to be how to let the code know which compositional field a point belongs to, and i don’t have much thoughts on this problem for this moment.

#4

First, I assume that you really meant to write this->introspection().name_for_compositional_index(i) instead of this->introspection().name_for_compositional_index(q) in your code snippet.

Second, yes, I think this looks at least conceptually correct. It is true that you will call name_for_compositional_index(i) many times this way and that that could be expensive. Whether it really is expensive is something you will only find out by profiling, but let’s at least assume for a moment that it is. Then I would suggest to just switch the order of the loops, i.e., use code along the lines of

  for(unsigned int i=0; i<this->introspection().n_compositional_fields(); ++i)
  {
    if(this->introspection().name_for_compositional_index(i) == “upper_crust”)
    {
      for (unsigned int q=0....)
        statement...
    }
    else if(this->introspection().name_for_compositional_index(i) == “lower_crust”)
    {
      for (unsigned int q=0....)
        statement...
    }
  ...

This way you pull the expensive computations out of the loop over quadrature points.

Does this make sense?
Best
W.