Coefficient &Q not initializing

[goinea]

Summary

I am creating a reference Q to send information to Eval in order to store in elmat_diff. I continue to get an error saying that I have used Q as a reference but not initialized it. Any hints?

Details

Code

elmat_diff = 0.; for (int j = 0; j < ir.GetNPoints(); j++) { // loop over quadrature points

const IntegrationPoint& ip = ir.IntPoint(j);
et->SetIntPoint(&ip);

// TODO: add the code for computation of element matrix corresponding to the diffusive term
fe->CalcDShape(ip, dshape);//added
et->SetIntPoint(&ip);//added
w = et->Weight();//added
w = ip.weight / (square ? w: w*w*w);//added
Mult(dshape, et->AdjugateJacobian(), dshapedxt);//added
w *= Q->Eval(et, ip);//added
AddMult_a_AAt(w, dshapedxt, elmat_diff);//added
//       and store it into "elmat_diff"

[mortezah]

Summary

I am creating a reference Q to send information to Eval in order to store in elmat_diff. I continue to get an error saying that I have used Q as a reference but not initialized it. Any hints?

Details

Code

elmat_diff = 0.; for (int j = 0; j < ir.GetNPoints(); j++) { // loop over quadrature points

const IntegrationPoint& ip = ir.IntPoint(j);
et->SetIntPoint(&ip);

// TODO: add the code for computation of element matrix corresponding to the diffusive term
fe->CalcDShape(ip, dshape);//added
et->SetIntPoint(&ip);//added
w = et->Weight();//added
w = ip.weight / (square ? w: w*w*w);//added
Mult(dshape, et->AdjugateJacobian(), dshapedxt);//added
w *= Q->Eval(et, ip);//added
AddMult_a_AAt(w, dshapedxt, elmat_diff);//added
//       and store it into "elmat_diff"

What is Q exactly?

Also, when creating an issue mention my name @mortezah . That way I will get a notification and will try to respond much quicker :)

[mortezah]

Also for the a4_advection_diffusion file I continue to get the singularity error: Apparently I have zeros in L, U, L+U, and LSUB, any hints. Code:

ifdef MFEM_USE_SUPERLU // if MFEM has super lu use a direct solve

SuperLUSolver superlu = new SuperLUSolver(MPI_COMM_WORLD); Operator SLU_A = new SuperLURowLocMatrix(A.As()); superlu->SetPrintStatistics(true); superlu->SetSymmetricPattern(false); superlu->SetColumnPermutation(superlu::METIS_AT_PLUS_A); superlu->SetRowPermutation(superlu::LargeDiag_MC64); superlu->SetIterativeRefine(superlu::SLU_DOUBLE); superlu->SetOperator(SLU_A); superlu->Mult(B, X);

I can think of 2 ways this could be caused.

1. Your shape functions or their derivatives are not implemented properly. You can rule this out if a4_diffusion gives you correct/reasonable results.
2. You are not setting the problem properly. Note that a4_advection_diffusion is very similar to a4_diffusion. So you can build things incrementally following exactly the implementation in a4_diffusion. For example, if turn off advective term (i.e. setting the vector a to zero) you should get the diffusion result. Do you get an error suggesting your matrix is singular without the advective term, too?

[goinea]

I have removed Q in seeing my error earlier, I do get a similar answer for diffusion. @mortezah

[mortezah]

I have removed Q in seeing my error earlier, I do get a similar answer for diffusion. @mortezah

Can you be more precise? Do you mean for diffusion problem you also get the same error? Or do you mean for diffusion case (by setting a = 0) you get the same answer as in a4_diffusion?

[goinea]

I get the same error. @mortezah

[mortezah]

I get the same error. @mortezah

do you get the same error when running a4_diffuison?

[goinea]

Yes @mortezah

[mortezah]

@goinea then, the most likely culprit is that your shape functions (or their derivatives) are not implemented properly.

Let's try the following:

1. when you run a4_projection it will create a vtk file. Can you take a screenshot of the field and paste it here?

[mortezah]

Attached is the screen shot

what I meant to say was to open that file in ParaView and then take the screenshot of the field visualized on the mesh.

[goinea]

@mortezah

[mortezah]

Okay, attached:

in the drop-down menu that says material pick the other option that says "field". That way you can look at the field value.

[mortezah]

@goinea that does not look correct. Note that you should see the magnitude of this https://github.com/SCOREC/fep/blob/d16f0b636671f519bcfb8ea91d07fc04638ad091/a4/a4_projection.cpp#L107 user-defined function. (you can also look at individual components by changing the value in the drop-down menu that currently says "Magnitude" in the picture you attached, to "X" or "Y".

I suggest going back to the implementations of your shape functions and making sure they are correct. There simple checks that you can do. For example, the shape function value has to be 1 at the node it corresponds to and zero at all other nodes.

[goinea]

Since I am in the bracketed section for say, Quads. Do I still have to denote the specific shape function: name: LG_QuadrilateralElement::CalcDShape...etc or can I just set CalcDShape(ip,etc)?

[goinea]

@mortezah

[mortezah]

you cannot call CalcShape or CalcDShape like that. You have to create a FiniteElement object first. Look at here https://github.com/SCOREC/fep/blob/d16f0b636671f519bcfb8ea91d07fc04638ad091/a4/a4_element_stiffness.cpp#L75-L77 for example. Once you have created that object (let's assume you are also calling it fe) then you can call CalcShape or CalcDShape using (note that I am not passing in all the arguments in the pseudo-codes below)

fe->CalcShape()
fe->CalcDShape()

[goinea]

I followed the form for fe->CalcShape(ip,etc) and my interpolation error is still not zero. Am I missing something or I just need to loop at some other literature?

[goinea]

@mortezah

[mortezah]

for order 1, what are the 4 shape function values that you get if ip is (0,0), (1,0), (0,1) and (1,1) i.e. the corners of a single reference quad.

[goinea]

1:(0,0),2:(1,0,)3:(1,1),4:(0,1) 1:0, 2:1/2, 3:1, 4:1/2

[goinea]

@mortezah

[mortezah]

1:(0,0),2:(1,0,)3:(1,1),4:(0,1) 1:0, 2:1/2, 3:1, 4:1/2

for each of those points, you will have 4 shape function values. I do not quite understand the numbers you provided.

[goinea]

@mortezah

[mortezah]

@goinea the notation that you are using seems to be for a quad defined in [-1,1]x[-1,1]. As it is stated here https://github.com/SCOREC/fep/blob/d16f0b636671f519bcfb8ea91d07fc04638ad091/a4/LagrangeElements.hpp#L198 in MFEM reference quads are defined on [0,1]x[0,1].

[goinea]

I thought I had adjusted for that, but how would that give me 4 values each? Best,

Sent from my iPhone

On Apr 15, 2021, at 8:42 PM, Morteza @.***> wrote:

 @goinea the notation that you are using seems to be for a quad defined in [-1,1]x[-1,1]. As it is stated here https://github.com/SCOREC/fep/blob/d16f0b636671f519bcfb8ea91d07fc04638ad091/a4/LagrangeElements.hpp#L198 in MFEM reference quads are defined on [0,1]x[0,1].

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[goinea]

@morteza

Sent from my iPhone

On Apr 15, 2021, at 8:42 PM, Morteza @.***> wrote:

 @goinea the notation that you are using seems to be for a quad defined in [-1,1]x[-1,1]. As it is stated here https://github.com/SCOREC/fep/blob/d16f0b636671f519bcfb8ea91d07fc04638ad091/a4/LagrangeElements.hpp#L198 in MFEM reference quads are defined on [0,1]x[0,1].

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[mortezah]

@goinea For order 1, you have 4 shape functions, one for each corner. Evaluating them at a given xi coordinate (ip.x, ip.y) should give you 4 values. At each corner you are expected to get one 1. and three 0.s (one would be for the shape function corresponding to that corner).

PS. You can mention @mortezah in the same comment to avoid having a separate comment with only @mortezah

[mortezah]

@goinea you are welcome.

As for references, I'm pretty sure most of these concepts are covered in class and references mentioned therein. You just have to be careful since different references might use different notations (example quads that are parametrized in different intervals). When designing a4, I tried to put lots of comments throughout the code to point out this potential differences.