FEAP User Forum
FEAP => Programming => Topic started by: a118145 on August 19, 2020, 06:46:22 AM
-
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:
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:
...
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
-
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?
-
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
-
Dear Professor Taylor,
I found the problem: In line 273 of the p_paraview routine, it says
write(plu,2000) sigl(ii)The corresponding format specifier in line 347 reads
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.
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
-
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.
-
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