Strange Projection Behaviour [v8.5]

[a118145]

Dear feap forum,

I recently encountered strange Paraview output files when using PVIEw,TIME. I implemented a linear elastic fictitious domain material routine, which I tested using a square domain with a circular inclusion. The attachment "orderParameter.png" shows the order parameter p with the circular inclusion p=0 and the surrounding material p=1. A smooth transition is established via a hypberbolic tangent and the gradient, which is needed for some computations in the vicinity of the interface is depicted by the white arrows. The order parameter and all gradient components are saved to the time independent history once at the beginning.

The input file reads as follows:

Code: [Select]

FEAP * * 2D: 6-node triangle
INCLude,PROB_2dellnco

MATErial,1
   user 01
     Umat 3
     E1 1000.0d0 ! stiff inclusion
     nu1 0.25d0     
     E2 100.0d0   ! matrix
     nu2 0.25d0

INCLude,MESH_2dellnco

EBOUndary
  1  0 1 0
  2  0 0 1
  1 10 1 0

EDISplacement
  1 10 2.0 0

END mesh

! Solution
BATCh
  PVIEw,TIME
  TOL,ITER,1.0d-7,1.0d-7
  LOOP iter 5
    TANG,,1
  NEXT iter
  PVIEw,TIME
END

STOP


As you can see, the paraview file is written twice and hence is the projection of Gauss point variables to the nodes. The first file (orderparameter.png) contains 12 values, which are zero (stress and strain), which is obvious, because loading hasn't started yet, and 4 nonzero values (order parameter and gradient). Everything looks fine at this moment. However, when I write the paraview file the second time, each of the 16 projected values shows strange behaviour as you can see from the exemplary picture I attached (strange_projection.png). To verify the simulation result, I had a look at the primary field (warped_displacement.png), but everything seems fine.

It gets even stranger, when I test the same setup with my linear elastic testing routine. For this case, I do not observe those oscillations/artifacts at all, but rather can observe the homogeneous stress and strain response. The relevant code snipped of the projection routine for the quadratic triangles reads as follows:

Code: [Select]

...
elseif (nel.eq.6 .and. ndm.eq.2) then ! quadratic tri elements
  mattri(:,:)  = 0.0d0
  siggtri(:,:) = 0.0d0
! loop over al QPs
  do l = 1,lint
    do i = 1,iste
      if (i.le.dep) then
        field(i)=hr(nh2+dep*(l-1)+i-1)
      else
        field(i)=hr(nh3+indep*(l-1)+(i-dep)-1)
      endif
    end do
    do j = 1,3
      xj          = el2(j,l)*el2(4,l)
      mattri(:,j)    = mattri(:,j)    + el2(1:3,l)*xj
      siggtri(j,1:iste) = siggtri(j,1:iste) + field(1:iste)*xj
    end do ! j
  end do ! l
  call invert(mattri,3,3)
  do j = 1,iste
    gsigtri(:,j) = mattri(:,1)*siggtri(1,j)+ mattri(:,2)*siggtri(2,j)+ mattri(:,3)*siggtri(3,j)
  end do ! j
  p(1:6) = 1.0d0

  ! Project stresses
  do j = 1,iste
    s(1:3,j) =  gsigtri(1:3,j)                   ! Vertex nodes
    s( 4,j)  = (gsigtri(1,j) + gsigtri(2,j))*0.5d0  ! Mid-edge node
    s( 5,j)  = (gsigtri(2,j) + gsigtri(3,j))*0.5d0
    s( 6,j)  = (gsigtri(3,j) + gsigtri(1,j))*0.5d0
  end do ! j
  write(*,*) dep,indep,iste,istv
...
The parameters dep and indep reflect the number of time dependent and independent history variables. The values iste and istv are set before the projection routine is called. For the linear elastic case, dep=12, indep=0, iste=istv=12 and for the fictitious domain case, dep=12, indep=4, iste=istv=16.

Does anybody have an idea, why the projection is sensitive to the number of projected values even if all variables seem to be set correctly?

Warm regards,

a118145

[Prof. R.L. Taylor]

What happens if you include 16 items in the first file too?  Is it possible paraview does not change dimensioning?

It is hard to say what is happening.  What do the plots in feap look like?  Do they also have the same type of behavior?

[a118145]

Dear Professor Taylor,

In all three files which I uploaded, 16 items are included. In the file orderParameter.png, the first 12 items are zero (which is expected) and 13 to 16 reflect the order parameter and its gradient.

I think, this is decrepit, because the history array does not change its dimension (see above):

Is it possible paraview does not change dimensioning?

Indeed I observed, that the plots in feap are correct. I suspect, that Paraview (or the PVIEw export routine) messes up the connectivity. The legend for the plots with artifacts pretty much covers the min and max values from the stresses (except for some very low values at around -300), they only seem "badly" matched to the nodes.

Warm regards,

a118145

[a118145]

Dear Professor Taylor,

I found the problem: In line 273 of the p_paraview routine, it says

Code: [Select]

write(plu,2000) sigl(ii)The corresponding format specifier in line 347 reads

Code: [Select]

2000  format(1p,1e13.5,$)
However, the specifier is not compatible to numbers, which have a very small exponent, such as 2.535E-320. What happens is, that the E is disappearing in the vtu file so that 2.535-320 is written to the vtu file. This format seems to mess Paraview up, because the -320 value,  which is approximately my min range of the legend, appears in all 16 entries. My fix for the problem is as follows. I modified line 267 such that not only the "none-zero"-ness is verified, but also the absolute value. If it is below 1e-16 (numerical zero), nothing is added, which keeps the vtu file clean.

Code: [Select]

if(abs(sig(jj,ii)).gt.1e-16) then
Surely, a better fix would be a more robust format specifier, but I leave this to you, because a short web search on this topic was not fruitful

Warm regards,

a118145

[Prof. R.L. Taylor]

This is another format issue where exponents are too large/small for default mode. 
You can also try setting the format as

2000  format(1p,1e14.5e3,$)

or if 5 digits are enough as

2000 format(1p,1e13.4e3,$)

The exponent -320 seems too small to be realistic, I am curious why?  Usually when I see this it comes from a program error of an integer interpreted as a real or similar.

[a118145]

Dear Professor Taylor,

this is indeed interesting. Good to know for the future.

The exponent -320 seems too small to be realistic, I am curious why?  Usually when I see this it comes from a program error of an integer interpreted as a real or similar.

I checked all user routines and spotted a typo in the strain variable. Thereby, some stress coordinates were set wrong and some others were not set at all. For the latter ones, the random allocation value X.YZe-320 was still present. The influence on the actual result was so small that you would not recognize it by looking at the paraview files. Now, everything is fine.

Warm regards and thank you again for your time,

a118145