Discussion Forum

Yu,

I would look into using the "Domain ODE's and DAE's" feature.

What this physics does is calculate (in a distributed fashion) a second order ODE. See the following blog post:

https://www.comsol.com/blogs/overview-integration-methods-space-time/

Scroll to the section "Temporal Integration by Means of Additional Physics Interfaces". You will have to define a new dependent variable that is solved for using the ODE interface, make its definition the variable you are wanting to integrate, and then put a 1 on the diagonal of the "Damping or Mass Coefficient" for that variable, and ensure the "Mass Coefficient" is 0 (since you are performing a first-order integration). This (when you compare to the equation listed on the interface) gives:

When the interface solves for "u", it will perform an integration from 0 to current time t of that variable and make it available to other physics interfaces by the dependent variable name you have defined.

If I haven't explained this clearly, perhaps someone else can chime in. Good luck!

Hi Yu,

As I understand it, you are asking for a spatial integration over each element, possibly for every time or parameter step. There is no direct way in which you could compute a result which is constant per element, and contains the integral of some expression over each element.

The main question here is the purpose of this. How you aregoing to use those results? They would for example be mesh size dependent. Maybe if you explain the purpose, it is possible to suggest a solution.

You can compute an approximation to what I think you are asking about without integration though. Use the expression 'centroid(u)*meshvol'. Rather than integrating the DOF 'u', you would then use the centroid value in each element and multiply it with the element volume. As long as the elements are not too distorted, and gradients through the element do not deviate much from linearity, the difference from true integration should be that large.

Computing the integral over a single specific element (say number 17) can be done using a standard integration operator with the expression 'u*(meshelement==17)'. Here the boolean expression will be zero in all but the target element.

Regards,
Henrik

Hi Yu,

As I understand it, you are asking for a spatial integration over each element, possibly for every time or parameter step. There is no direct way in which you could compute a result which is constant per element, and contains the integral of some expression over each element.

The main question here is the purpose of this. How you aregoing to use those results? They would for example be mesh size dependent. Maybe if you explain the purpose, it is possible to suggest a solution.

You can compute an approximation to what I think you are asking about without integration though. Use the expression 'centroid(u)*meshvol'. Rather than integrating the DOF 'u', you would then use the centroid value in each element and multiply it with the element volume. As long as the elements are not too distorted, and gradients through the element do not deviate much from linearity, the difference from true integration should be that large.

Computing the integral over a single specific element (say number 17) can be done using a standard integration operator with the expression 'u*(meshelement==17)'. Here the boolean expression will be zero in all but the target element.

Regards,
Henrik

Hi, Henrik:

Thank you so much for your reply. Your answer is exactly what I want to get. I want to calculate a volume averaged value for a variable on element level, and then use it in a weak form. To my understand, the centroid operator you introduced here should be the volume averaged value, say (f u dVe)/Ve, here f is the integrate symbol. I am not sure whether I am correct. Would you mind further explaining the centroid operator at a more mathematical level?

Thanks a lot for your help.

Yu

Hi Yu,

The centroid() operator just picks the value at the centroid of the element, so it is not a true average. But as long as the solution is reasonably smooth, it will give a decent approximation.

An even better approximation can be found using the expression 'gpeval(4,u,0)'. This expression picks the value of 'u' in the (fourth order) Gauss points of the element, and makes a least-squares fit to a constant value.

If you are interested in assesing the accuracy of these approaches, you can for example plot 'centroid(u)-gpeval(4,u,0)' or compare the integrals over the whole domain of 'u', 'centroid(u)', and 'gpeval(4,u,0)' respectively.

You can read more about these operators in the user's guide.

Regards,
Henrik

Hi Yu,

One of my colleagues came up with a much better solution (thanks Dan!).

This requires an integration operator created under Definitions, say intop1().

Then ‘intop1((meshelement==dest(meshelement))*u’ will give the true (as far as the numerics goes) per element integral of 'u'.

To compute the true per element average, use ‘intop1((meshelement==dest(meshelement))*u)/intop1((meshelement==dest(meshelement)))’

Regards,
Henrik

Hi Yu,

One of my colleagues came up with a much better solution (thanks Dan!).

This requires an integration operator created under Definitions, say intop1().

Then ‘intop1((meshelement==dest(meshelement))*u’ will give the true (as far as the numerics goes) per element integral of 'u'.

To compute the true per element average, use ‘intop1((meshelement==dest(meshelement))*u)/intop1((meshelement==dest(meshelement)))’

Regards,
Henrik

Hi, Henrik:

That's fantastic! Now I can see how powerful COMSOL is. Thanks a lot for your such a detailed and valuable reply.

Best regards,

Yu

Hi Yu,

One of my colleagues came up with a much better solution (thanks Dan!).

This requires an integration operator created under Definitions, say intop1().

Then ‘intop1((meshelement==dest(meshelement))*u’ will give the true (as far as the numerics goes) per element integral of 'u'.

To compute the true per element average, use ‘intop1((meshelement==dest(meshelement))*u)/intop1((meshelement==dest(meshelement)))’

Regards,
Henrik

Hi, Henrik:

Thank you again for your attention to my message. I have tried the three approaches you provided with me. I found that the approach "intop1((meshelement==dest(meshelement))*u" is exremely computionally expansive when being used in weak form. The operator "centroid" gives a more efficient way, but as you mentioned it is just a value on the centroid. The operator "gpeval" gives a decent result, but still takes much more computational time. My objective of introducing the volume average term in weak form is to increase the computional efficiency of the model. So I am wondering whether there is another more robust and efficient way to determine the volume average value? Is it possible to do that by programming through Matlab?

Thank you so much for your help and time.

Best regards,

Yu

Hi Yu,

Would not using original non-averaged field be both the most accurate and the most efficient approach?

Why does your weak expression need the average? Is the weak expression coupling to a field of lower order? Do you want to reduce the number integration points for the weak contribution?

By the way, you must consider wether this field is just a source in your weak expression, or wheter the coupling works both ways.

Regards,
Henrik

Hi, Henrik:

The introduction of average value in the weak form is to create a stabilized technique, like Isotropic diffusion technique in COMSOL to the stabilize the convection-dominated transport problems. It plays a role of artifical diffustion (source term) for the whole model.

Are there any other possible ways to get the element average value? I think the centroid value in COMSOL should be determined from one-point Gauss quadrature at the centroid point. So, it can be regarded as a reduced Gauss quadrature, right? I am nore sure whether I am correct. If correct, I think the centroid value may be able to represent the volume average value to some extent.

Thank you for your time and help.

Yu

Hi Yu,

Using the centroid value is similar, but not necessarily equivalent, to using a 1 point integration.

If the expression you integrate is of the type centroid(u)*test(v), then test(v) is only evaluated at the centroid if you use 1 point integration. Also, this means that for example the element volume is represented by its scale (Jacobian) only at the centroid.

As for the methods for evaluting the element average, I think all options have been covered above.

Regards,
Henrik