Topic: global2local

Hi,

In a relatively recent commit (6143c82ef3c8f0f8ce5c2eda01e9b2c5b516a37a "Fixing global2local implementation"),
fei2dtrquad.C, fei2dquadquad.C, fei3dtetquad.C och fei3dwedgelin.C, the global2local function was changed
from jac.solveForRhs(res, delta, true); to jac.solveForRhs(res, delta);
This change removes a transpose from the Jacobian.

This change destroys the convergence in global2local for a problem I am currently working with.
My question is therefore:
Has anyone verified that it was indeed correct to modify the transpose in the functions mentioned above?

Performing the derivation myself I get the Jacobian as
J_ab = sum_i x_a^i (dN^i/dxi)_b,

but in  giveJacobianMatrixAt, the Jacobian is computed as
J_ab =  (dN^i/dxi)_a sum_i x_b^i.

Hence, it seems to me that the transpose should actually be there.

Regards,
Erik

2 (edited by nitramkaroh 01-11-2016 19:08:54)

Re: global2local

Dear Erik,

Some time ago, I found a problem in the function global2local in the class fei3dhexaquad.C and it was corrected by removing the transposition.
Recently, I had the same problem with global2local in fei3dwedgelin.C, so I changed it in this file as well as in other classes. I checked the files again, and it seems to me that the problem is caused by the inconsistent definition of the Jacobian. If you look into fei3dwedgelin or fei3dhexaquad, it seem to me that the Jacobian is computed as J_ab = sum_i x_a^i (dN^i/dxi)_b, while in fei2dquadquad as J_ab =  (dN^i/dxi)_a sum_i x_b^i.

Martin

Re: global2local

Hi Martin,

Thank you for the quick reply.  You are right, it seems like the Jacobian is computed differently for 2D and 3D elements.
I think the ideal solution would be to change to J_ab = sum_i x_a^i (dN^i/dxi)_b in the 2D elements as well.
However, giveJacobianMatrixAt is called by some other functions as well, so simply changing in giveJacobianMatrixAt could affect other parts of the code, and I don't know if these functions are all tested. I can try to take a look at this.

Regards,
Erik

Re: global2local

Going through the uses (did I miss something?)

oofemlib/meshqualityerrorestimator.C:94:        fei.giveJacobianMatrixAt( jac, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(elem) );
doesn't care about the transpose.

sm/Elements/3D/lspace.C:391:     * this->interpolation.giveJacobianMatrixAt (jm, domain, nodeArray, coords);
code is actually commented, but it looks like it was trying to use jm as a way to define the local c.s.

sm/Elements/structural3delement.C:130:        this->giveInterpolation()->giveJacobianMatrixAt( jac, lcoords, FEIElementGeometryWrapper(this) );
sm/Elements/structural2delement.C:156:        this->giveInterpolation()->giveJacobianMatrixAt( jac, lcoords, * this->giveCellGeometryWrapper() );
one of these are wrong. I wrote these and I basically try to extract the tangent (for defining material orientation).
Here I assumed

J = [ dx/dxi, dy/dxi
      dx/deta, dy/deta]

sm/Elements/Shells/quad1mindlinshell3d.C:269:    this->interp.giveJacobianMatrixAt(jac, localCoords, FEIVertexListGeometryWrapper(lnodes) );
this was untested experimental code. It would have to be tested anyway, so you can ignore this.

sm/Elements/Shells/solidshell.C:309:    interp->giveJacobianMatrixAt(J0mat, center, FEIElementGeometryWrapper(this) );
not sure about this one. It does pick apart the the columns and defining a coordinate system using them, but I can't quite tell what it actually wants.
Jim wrote this, so you could ask him directly Erik, or we assume that the current way it does it is what is truly wants (it expects the 3D version to be true).

Re: global2local

I made a default implementation for global2local in FEInterpolation2d and pushed it to Mickes repo.
Have a look and tell me what you think.
/Erik

Re: global2local

There are actually some small differences in triangle elements, , i.e. should the output be 2 or 3 values.
Now, it'll be 2, but I think it used to be 3. Of course, the third value is redundant, and we should consider dropping it universally.


I'm also quite sure that
FEI2dTrLin :: inside
is wrong, calculations for triangles should be very different from quads.

and

+    if ( error > convergence_limit ) { // Imperfect, could give false negatives.
+        OOFEM_WARNING("Failed convergence");
+        answer = {1. / 3., 1. / 3.};
+        return false;
+    }

is not correct for quads.  I don't even know if we should return anything here. Maybe just clear it instead, since something has gone wrong.

Re: global2local

Yes, checking whether the point is inside is different for triangles and quads. Therefore, I made an inside() function that is virtual and overloaded by FEI2dTrLin, FEI2dQuadLin, ....
I do believe that the inside check is correct, have I missed something here?

Regarding the second point: I agree, since something has gone wrong, we could clear it, or even throw an error.

Re: global2local

Yes, i think you missed something.

bool FEI2dTrQuad :: inside(const FloatArray &lcoords) const
{
    const double point_tol = 1.0e-3;
    bool inside = true;
    for ( int i = 1; i <= 2; i++ ) {
        if ( lcoords.at(i) < ( -1. - point_tol ) ) {
            inside = false;
        } else if ( lcoords.at(i) > ( 1. + point_tol ) ) {
            inside = false;
        }
    }

It does the same check as for quads, which is not correct. It should check 0 <= lcoords <= 1 and lcoords(0) + lcoords(1) <= 1

Re: global2local

You are right. I have now corrected the bounds in the inside check.