Nodal projection of user variable calculated in integration points

[JReinoso]

Dear All,

here is the problem. I am attaching a couple of files:

- Icohe3D, which is an example of 2 hexahedrical elements with an interface in the model. The solid elements are FEAP isotropic elastic elements whereas the interface model is a cohesive user element.
- elmt47.f is the cohesive user element.

My idea is to plot the cohesive damage variable from the user element. As it is computed in the integration points, like the stresses, I have thought of doing the nodal projection (following slcn2d.f) of this magnitude and then use the paraview file to plot it (following http://feap.berkeley.edu/forum/index.php?topic=1960.0).

You can see in elmt47.f that the projection is done in the gaussValueProjection subroutine, where I define the damage variable in the slot 16 (as FEAP in solid 3D elements creates 15 by default). Subsequently, when I run the file, I check on Paraview that the slot is created but somehow is not being updated with the cohesive damage magnitude. I wonder if this is due to combining a FEAP element with an user element.

Any help on this is much appreciated. Many thanks.
JReinoso.

[Prof. S. Govindjee]

Instead of

Code: [Select]

STREss,ALL try

Code: [Select]

STREss,NODE,1 I think you need to force the projection to the nodes, something that STREss,ALL does not do.  STREss,NODE will force the projection; alternately if you plot any stresses using FEAP's builti-in plotting it will force the projection.

[JReinoso]

Dear Prof. Govindjee,

Thanks for your quick reply. I tried changing that line you mentioned in the example and unfortunately, I am still getting a blank isocontour. I can tell you that in the gaussValueprojection subroutine, the p(nel) and s(nel,*) matrices are being updated accordingly with the cohesive damage but somehow I do not manage to translate them into the postprocessing (using both built-in FEAP and Paraview visualization tools).

JReinoso.

[FEAP_Admin]

A couple of comments.

(1) Under isw.eq.1 when you set istv you should do it as

Code: [Select]

istve=max(16,istv)in case there are other elements in your mesh that need more slots.

(2) In your gaussValueProjection routine, it looks like you are putting stress-1 into plotting stress-16.  You probably want something more like

Code: [Select]

xg = xsj(l)*shp3(i,l)
r(i) = r(i) + xg
s(i,1:16) = s(i,1:16) + sig(1:16,l)*xg

[Prof. R.L. Taylor]

It looks like you are passing "nes" as "1" when you call gaussValueProjection.  If you only want to project the one value you need to ensure that  all the numbers match correctly since the stress is dimensioned sig(1,*) with what you have. What is the real dimension for the stress and other arrays??

[JReinoso]

It looks like you are passing "nes" as "1" when you call gaussValueProjection.  If you only want to project the one value you need to ensure that  all the numbers match correctly since the stress is dimensioned sig(1,*) with what you have. What is the real dimension for the stress and other arrays??

Dear Prof. Taylor,

Thanks for your reply. Yes indeed, nes is 1, as I only want to plot the cohesive damage. Therefore, sig should be (1,*). As I want to plot it in the 16th slot, I do s(1:nel,16). However, since the other 15 are left blank, I have modified the code to plot the sig(1,*) into s(1:nel,1), the first slot. Still, I cannot manage to plot the cohesive damage neither in the built-in FEAP visualization tool nor in Paraview.

A couple of comments.

(1) Under isw.eq.1 when you set istv you should do it as

Code: [Select]

istve=max(16,istv)in case there are other elements in your mesh that need more slots.

(2) In your gaussValueProjection routine, it looks like you are putting stress-1 into plotting stress-16.  You probably want something more like

Code: [Select]

xg = xsj(l)*shp3(i,l)
r(i) = r(i) + xg
s(i,1:16) = s(i,1:16) + sig(1:16,l)*xg

Dear FEAP_Admin,

Thanks for your reply. I have tried doing both things (as I mentioned Prof. Taylor, I have modified the code to have only one slot. As nes = 1, I try to plot it in s(1:nel,1)), but still I cannot manage to have the results plotted.

In the new code, I am writing p(1:nel) and s(1:nel,*) into the terminal and as you can see, the results are being properly added into the matrices. Therefore, my guess is that the problem has to be in the input file.

JReinoso.

[Prof. R.L. Taylor]

When I look at elmt47.f and what happens when isw.eq.8 there does not seem to be anywhere where the "sig" is computed.  There must be a value for each of the quadrature points you are projecting.  Look at how, say sld2d1.f, does it (./elements/solid2d/sld2d1.f).  If you have the values somewhere then use them (I see something called "sigforc" in and older call, is that what you want???).  I would suggest printing everything to the output file (iow) to see what is missing.

[JReinoso]

Dear All,

The problem is now solved. I was missing the definition of the pstyp variable in the isw.eq.1 and now it's working properly. Thanks a lot for the help provided, though!

Best regards,
JReinoso.