So, I tried again with this problem to see if I can solve it using the provided input file. The hint regarding printing out the stresses actually helped a lot.
STRE, ELEM for the ntie elements output for Ibeam (the sample input file):
Feap * * L-shaped cantilever beam
Interface Element 27 Results
Elmt 1-Coord 2-Coord 3-Coord Lambda Constraint Angle
51 5.000E+00 1.000E+01 0.000E+00 0.0000E+00 -2.2314E-07 9.0000E+01
51 5.000E+00 1.000E+01 3.333E+00 -3.1070E-03 4.2414E-07 9.0000E+01
52 5.000E+00 1.000E+01 3.333E+00 -3.1070E-03 -2.1207E-07 9.0000E+01
52 5.000E+00 1.000E+01 1.000E+01 -3.4582E-04 7.5083E-07 9.0000E+01
53 5.000E+00 1.000E+01 1.000E+01 -3.4582E-04 -5.0055E-07 9.0000E+01
53 5.000E+00 1.000E+01 2.000E+01 -3.5300E-03 4.1217E-07 9.0000E+01
54 5.000E+00 1.000E+01 2.000E+01 -3.5300E-03 -4.1217E-07 9.0000E+01
54 5.000E+00 1.000E+01 3.000E+01 -3.0809E-02 2.6722E-07 9.0000E+01
55 5.000E+00 1.000E+01 3.000E+01 -3.0809E-02 -2.6722E-07 9.0000E+01
55 5.000E+00 1.000E+01 4.000E+01 -1.4404E-01 -2.0457E-08 9.0000E+01
56 5.000E+00 1.000E+01 4.000E+01 -1.4404E-01 3.0685E-08 9.0000E+01
56 5.000E+00 1.000E+01 4.667E+01 -5.6433E-01 -3.7679E-08 9.0000E+01
57 5.000E+00 1.000E+01 4.667E+01 -5.6433E-01 7.5358E-08 9.0000E+01
57 5.000E+00 1.000E+01 5.000E+01 -1.7135E+00 -6.9918E-16 9.0000E+01
which seems to work properly. For my input file, at first, I got:
Interface Element 27 Results
Elmt 1-Coord 2-Coord 3-Coord Lambda Constraint Angle
65 5.000E-01 0.000E+00 5.000E-01 -5.3401E-10 0.0000E+00 9.0000E+01
65 5.000E-01 6.250E-02 4.375E-01 3.7953E-01 9.7963E-01 1.1585E+01
66 5.000E-01 6.250E-02 5.000E-01 3.7953E-01 0.0000E+00 9.0000E+01
66 5.000E-01 1.875E-01 4.375E-01 4.9367E-01 9.8735E-01 9.1237E+00
67 5.000E-01 1.875E-01 5.000E-01 4.9367E-01 0.0000E+00 9.0000E+01
67 5.000E-01 3.125E-01 4.375E-01 4.9386E-01 9.8772E-01 8.9867E+00
68 5.000E-01 3.125E-01 5.000E-01 4.9386E-01 0.0000E+00 9.0000E+01
68 5.000E-01 4.375E-01 4.375E-01 4.9398E-01 9.8796E-01 8.9002E+00
69 5.000E-01 4.375E-01 5.000E-01 4.9398E-01 0.0000E+00 9.0000E+01
69 5.000E-01 5.625E-01 4.375E-01 4.9405E-01 9.8810E-01 8.8473E+00
70 5.000E-01 5.625E-01 5.000E-01 4.9405E-01 0.0000E+00 9.0000E+01
70 5.000E-01 6.875E-01 4.375E-01 4.9411E-01 9.8821E-01 8.8058E+00
71 5.000E-01 6.875E-01 5.000E-01 4.9411E-01 0.0000E+00 9.0000E+01
71 5.000E-01 8.125E-01 4.375E-01 4.9418E-01 9.8835E-01 8.7541E+00
72 5.000E-01 8.125E-01 5.000E-01 4.9418E-01 0.0000E+00 9.0000E+01
72 5.000E-01 9.375E-01 4.375E-01 6.0557E-01 9.8857E-01 8.6719E+00
73 5.000E-01 9.375E-01 5.000E-01 6.0557E-01 0.0000E+00 9.0000E+01
73 5.000E-01 1.000E+00 4.375E-01 9.8240E-01 9.8240E-01 1.0765E+01
Clearly, something was not correct as my problem is a 2D one, and should not have non-zero z-coordinates. So, I had a closer look into the code and it seems like some parts of stif27.f and stre27.f, where coordinates are used/calculated, are hard-coded with ndm=3. For example,
Ta(1:3,1) = 0.0d0
Ta(1:3,2) = 0.0d0
do i = 1,nel1
Ta(1:3,1) = Ta(1:3,1) + xa(1:3,i)*shpa(1,i)
Ta(1:3,2) = Ta(1:3,2) + xa(1:3,i)*shpa(2,i)
end do ! i
where although xa is dimensioned xa(ndm,4), the range for the first dimension of xa is not 1:ndm.
I modified the code to make it consistent with ndm. This time I get:
Interface Element 27 Results
Elmt 1-Coord 2-Coord 3-Coord Lambda Constraint Angle
65 5.000E-01 0.000E+00 0.000E+00 0.0000E+00 0.0000E+00 9.0000E+01
65 5.000E-01 6.250E-02 0.000E+00 2.3934E-12 0.0000E+00 9.0000E+01
66 5.000E-01 6.250E-02 0.000E+00 2.3934E-12 0.0000E+00 9.0000E+01
66 5.000E-01 1.875E-01 0.000E+00 9.3823E-12 0.0000E+00 9.0000E+01
67 5.000E-01 1.875E-01 0.000E+00 9.3823E-12 0.0000E+00 9.0000E+01
67 5.000E-01 3.125E-01 0.000E+00 9.4461E-12 0.0000E+00 9.0000E+01
68 5.000E-01 3.125E-01 0.000E+00 9.4461E-12 0.0000E+00 9.0000E+01
68 5.000E-01 4.375E-01 0.000E+00 9.9765E-12 0.0000E+00 9.0000E+01
69 5.000E-01 4.375E-01 0.000E+00 9.9765E-12 0.0000E+00 9.0000E+01
69 5.000E-01 5.625E-01 0.000E+00 1.0069E-11 0.0000E+00 9.0000E+01
70 5.000E-01 5.625E-01 0.000E+00 1.0069E-11 0.0000E+00 9.0000E+01
70 5.000E-01 6.875E-01 0.000E+00 1.0831E-11 0.0000E+00 9.0000E+01
71 5.000E-01 6.875E-01 0.000E+00 1.0831E-11 0.0000E+00 9.0000E+01
71 5.000E-01 8.125E-01 0.000E+00 1.0273E-11 0.0000E+00 9.0000E+01
72 5.000E-01 8.125E-01 0.000E+00 1.0273E-11 0.0000E+00 9.0000E+01
72 5.000E-01 9.375E-01 0.000E+00 2.2884E-11 0.0000E+00 9.0000E+01
73 5.000E-01 9.375E-01 0.000E+00 2.2884E-11 0.0000E+00 9.0000E+01
73 5.000E-01 1.000E+00 0.000E+00 2.4000E-11 0.0000E+00 9.0000E+01
But, I think these results are also not correct as this time the "constraint" value is zero, which seems not right. Also, the overall results of 2-patches are still different from the 1-patch setting.
I also tried changing ndm for my input file to ndm=3 to see if it works or not, but somehow I get zero (global) derivatives of the shape function in my user-element (but, this, I guess, is a different issue).
So, the open questions would be:
- Are the changes that I have done in the code sufficient to make it possible to use the slope element for ndm = 2? (a complete list of changes is provided below)
- Or, should I try to fix the problem of my user element with ndm=3, and use the original formulation for the slope element?
--- a/ver85/igafeap/elements/slope/stif27.f
+++ b/ver85/igafeap/elements/slope/stif27.f
@@ -180,2 +180,2 @@
- Ta(1:3,1) = Ta(1:3,1) + xa(1:3,i)*shpa(1,i)
- Ta(1:3,2) = Ta(1:3,2) + xa(1:3,i)*shpa(2,i)
+ Ta(1:ndm,1) = Ta(1:ndm,1) + xa(1:ndm,i)*shpa(1,i)
+ Ta(1:ndm,2) = Ta(1:ndm,2) + xa(1:ndm,i)*shpa(2,i)
@@ -187,2 +187,2 @@
- Tb(1:3,1) = Tb(1:3,1) + xb(1:3,i)*shpb(1,i)
- Tb(1:3,2) = Tb(1:3,2) + xb(1:3,i)*shpb(2,i)
+ Tb(1:ndm,1) = Tb(1:ndm,1) + xb(1:ndm,i)*shpb(1,i)
+ Tb(1:ndm,2) = Tb(1:ndm,2) + xb(1:ndm,i)*shpb(2,i)
@@ -238,2 +238,2 @@
- Ta(1:3,1) = Ta(1:3,1) + xsa(1:3,i)*shpa(1,i)
- Ta(1:3,2) = Ta(1:3,2) + xsa(1:3,i)*shpa(2,i)
+ Ta(1:ndm,1) = Ta(1:ndm,1) + xsa(1:ndm,i)*shpa(1,i)
+ Ta(1:ndm,2) = Ta(1:ndm,2) + xsa(1:ndm,i)*shpa(2,i)
@@ -244,2 +244,2 @@
- Tb(1:3,1) = Tb(1:3,1) + xsb(1:3,i)*shpb(1,i)
- Tb(1:3,2) = Tb(1:3,2) + xsb(1:3,i)*shpb(2,i)
+ Tb(1:ndm,1) = Tb(1:ndm,1) + xsb(1:ndm,i)*shpb(1,i)
+ Tb(1:ndm,2) = Tb(1:ndm,2) + xsb(1:ndm,i)*shpb(2,i)
@@ -431,4 +431,4 @@
- c_Ta(1:3,1) = c_Ta(1:3,1)
- & + (xa(1:3,i)+uca(1:3,i))*shpa(1,i)
- c_Ta(1:3,2) = c_Ta(1:3,2)
- & + (xa(1:3,i)+uca(1:3,i))*shpa(2,i)
+ c_Ta(1:ndm,1) = c_Ta(1:ndm,1)
+ & + (xa(1:ndm,i)+uca(1:ndm,i))*shpa(1,i)
+ c_Ta(1:ndm,2) = c_Ta(1:ndm,2)
+ & + (xa(1:ndm,i)+uca(1:ndm,i))*shpa(2,i)
@@ -439,4 +439,4 @@
- c_Tb(1:3,1) = c_Tb(1:3,1)
- & + (xb(1:3,i)+ucb(1:3,i))*shpb(1,i)
- c_Tb(1:3,2) = c_Tb(1:3,2)
- & + (xb(1:3,i)+ucb(1:3,i))*shpb(2,i)
+ c_Tb(1:ndm,1) = c_Tb(1:ndm,1)
+ & + (xb(1:ndm,i)+ucb(1:ndm,i))*shpb(1,i)
+ c_Tb(1:ndm,2) = c_Tb(1:3,2)
+ & + (xb(1:ndm,i)+ucb(1:ndm,i))*shpb(2,i)
diff --git a/ver85/igafeap/elements/slope/stre27.f b/ver85/igafeap/elements/slope/stre27.f
index 60e4015..99d4d5a 100644
--- a/ver85/igafeap/elements/slope/stre27.f
+++ b/ver85/igafeap/elements/slope/stre27.f
@@ -156,2 +156,2 @@
- Ta(1:3,1) = Ta(1:3,1) + xa(1:3,i)*shpa(1,i)
- Ta(1:3,2) = Ta(1:3,2) + xa(1:3,i)*shpa(2,i)
+ Ta(1:ndm,1) = Ta(1:3,1) + xa(1:ndm,i)*shpa(1,i)
+ Ta(1:ndm,2) = Ta(1:3,2) + xa(1:ndm,i)*shpa(2,i)
@@ -163,2 +163,2 @@
- Tb(1:3,1) = Tb(1:3,1) + xb(1:3,i)*shpb(1,i)
- Tb(1:3,2) = Tb(1:3,2) + xb(1:3,i)*shpb(2,i)
+ Tb(1:ndm,1) = Tb(1:ndm,1) + xb(1:ndm,i)*shpb(1,i)
+ Tb(1:ndm,2) = Tb(1:ndm,2) + xb(1:ndm,i)*shpb(2,i)
@@ -214,2 +214,2 @@
- Ta(1:3,1) = Ta(1:3,1) + xsa(1:3,i)*shpa(1,i)
- Ta(1:3,2) = Ta(1:3,2) + xsa(1:3,i)*shpa(2,i)
+ Ta(1:ndm,1) = Ta(1:ndm,1) + xsa(1:ndm,i)*shpa(1,i)
+ Ta(1:ndm,2) = Ta(1:ndm,2) + xsa(1:ndm,i)*shpa(2,i)
@@ -220,2 +220,2 @@
- Tb(1:3,1) = Tb(1:3,1) + xsb(1:3,i)*shpb(1,i)
- Tb(1:3,2) = Tb(1:3,2) + xsb(1:3,i)*shpb(2,i)
+ Tb(1:ndm,1) = Tb(1:ndm,1) + xsb(1:ndm,i)*shpb(1,i)
+ Tb(1:ndm,2) = Tb(1:ndm,2) + xsb(1:ndm,i)*shpb(2,i)
@@ -268 +268 @@
- xx(:) = xx(:) + xa(1:3,ixl(a))*shm(a)
+ xx(1:ndm) = xx(1:ndm) + xa(1:ndm,ixl(a))*shm(a)
(END)