Problem with the bending strip element

[resammc]

Dear all,

I need a multi-patch two-dimensional domain with C1-continuity across single patches for my calculations. The first thing which came to my mind was to try the Bending Strip element implemented in igaFEAP. It seems like it can be used for my purpose. However, I cannot make it to work properly.

To test the implementation, I use a square model and compare the stresses for a C1-continuous single patch domain versus a multi-patch domain consisting of 2 patches where the bending strip is working on the C0-continuous kink between the patches. The studies are done on a fairly coarse bi-quadratic NURBS mesh for a plane strain setting (this is to be able to validate the results using FEAP standard elements).

Since the strip works as a series of 6-node elements, I define an extra patch for it where the parametric direction transverse to the intersection is quadratic (only one element in this direction), and the other direction is linear (which covers the intersection between the patches). The control points are coming from the original mesh.

The problem is that, even for this simple case, the stresses are not distributed in the same way in the 2-patches case with the bending strip as the stresses of the 1-patch case.

So, here are my questions:

- I think I'm not using the nurbs_slope element correctly. Could you please have a look and let me know if something is not correct in the input file?

- Is it possible to directly use this element for planar problems like this?

- Is the orientation of the patch important? (in other words, is it important if the quadratic direction is the first parametric direction or the second?)

- Is the implementation based on the "The bending strip method for isogeometric analysis of Kirchhoff–Love shell structures comprised of multiple patches" paper by Kiendl et al.? Or is there another source that has explained it in more detail? In particular, I don't understand "ie1  / 1,2,4,3 /" and "ie2  / 2,1,5,6 /" (seems like it is not a usual control point numbering)

Kindly find attached my input files (the boundary conditions are chosen so that it is possible to have a non-uniform stress distribution to see the differences between the models).

Any help in this regard is greatly appreciated.

[Prof. R.L. Taylor]

Have you tried the NTIE command?  The patches need to match to use. See Sect 2.6.1 of IgA manual

[resammc]

Dear Professor Taylor,

thank you very much for your response, and for your hint. Somehow, I overlooked that section in the manual... it makes using the feature much easier and answers my question regarding the node numbering. And yes, the patches were compatible.

However, I still cannot get the bending strip to work correctly. It probably is easier to see the difference with my user element (a strain gradient model), so the figures are attached.

I use:

Code: [Select]

mate 2
NURB_SLOPE
QUADRATURE GAUSS
PENALTY,,penalty_value
for the material card. For the penalty value I tried 1.0d-1 to 1.0d-6 (as it seems like the element uses the inverse of this value), but I don't see any noticeable changes. I could not use the LAGRANGE multipliers option at all (no convergence).

[Prof. R.L. Taylor]

Solution uses perturbed Lagrangian. So penalty is large

[resammc]

I tried various values, but unfortunately, I cannot make it right. Is there any sample input file against which I can compare my input file?

[Prof. R.L. Taylor]

Attached is a simple example.

Also, if you output the "stre" for the "ntie" elements and post it may give a clue to what is happening.

[resammc]

Thank you very much for the hint and the input file. Unfortunately, I don't have access to my test system at the moment. I will try the input file and have a look at the stress outputs as you suggested and will report back here.

[resammc]

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):

Code: [Select]

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:

Code: [Select]

     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,

Code: [Select]

        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:

Code: [Select]

     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?

Code: [Select]

--- 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)


[Prof. R.L. Taylor]

What input file are you using?  Can you post it and I will try to debug.  There are a few things I see that need to be checked better.

[resammc]

Kindly find attached the correspnding input files.

[Prof. R.L. Taylor]

The current bending strip is not intended for a 2-d problem.  It is deduced for shells, and then replaced by the ntie command.

In a 2-d problem what you need to do is to have a form that ensures that the continuity in the normal direction to the edge satisfies  the continuity of the gradient.  So the condition is that du/dn are equal.  For equal meshing is is probably fairly easy but for a general case where control points lie on a straight line it involves the knot spacing in the derivative.

It appears that for the attached figure the increment in displacements u_62 - u_6 = u_6 - u_5  -> u_62 - 2*u_6 + u_5 = 0 is the correct constraint (it is the finite difference constant curvature).  Etc. for others.

[Prof. R.L. Taylor]

Figure

[resammc]

Thank you very much for the hints. I'll try to make it work...