Fixes in magnetic BC

Removed part of yesterday's changes to interpolation classes, which turned out to be unnecessary as of today

src/oofemlib/fei3dhexalin.C

@@ -572,53 +572,6 @@ FEI3dHexaLin :: surfaceLocal2global(FloatArray &answer, int iedge,
                    n.at(3) * cellgeo.giveVertexCoordinates( nodes.at(3) ).at(3) + n.at(4) * cellgeo.giveVertexCoordinates( nodes.at(4) ).at(3);
 }
 
-void FEI3dHexaLin::surfaceLocal2fullLocal( FloatArray &answer, int iSurf, const FloatArray &surfacelcoords, const FEICellGeometry &cellgeo ) const
-//based on iSurf, calculates the full local coordinate array
-{
-    answer.resize( 3 );
-
-    //the ordering and flipping of the local coordinates has been determined empirically on simple cube example
-    //by comparison of global coordinates returned from surface and full coordinate arrays
-    switch ( iSurf ) {
-    case 1:
-        // zeta = 1
-        answer = { -surfacelcoords.at( 1 ), -surfacelcoords.at( 2 ), 1.0 };
-        break;
-    case 2:
-        // zeta = -1
-        answer = { -surfacelcoords.at( 2 ), -surfacelcoords.at( 1 ), -1.0 };
-        break;
-    case 3:
-        // ksi = -1
-        answer = { -1.0, -surfacelcoords.at( 1 ), surfacelcoords.at( 2 ) };
-        break;
-    case 4:
-        // eta = 1
-        answer = { -surfacelcoords.at( 1 ), 1.0, surfacelcoords.at( 2 ) };
-        break;
-    case 5:
-        // ksi = 1
-        answer = { 1.0, surfacelcoords.at( 1 ), surfacelcoords.at( 2 ) };
-        break;
-    case 6:
-        // eta = -1
-        answer = { surfacelcoords.at( 1 ), -1.0, surfacelcoords.at( 2 ) };
-        break;
-    }
-
-    //debug: check whether this works
-   /* FloatArray globalFromSurface, globalFromFull, error;
-    this->surfaceLocal2global( globalFromSurface, iSurf, surfacelcoords, cellgeo );
-    this->local2global( globalFromFull, answer, cellgeo );
-
-    error = globalFromSurface - globalFromFull;
-
-    OOFEM_LOG_INFO("\n\nIsurf %i -- surface/global/error\n --", iSurf);
-    globalFromSurface.printYourself();
-    globalFromFull.printYourself();
-    error.printYourself();*/
-}
-
 double
 FEI3dHexaLin :: surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
 {

src/oofemlib/fei3dhexalin.h

@@ -81,7 +81,6 @@ class OOFEM_EXPORT FEI3dHexaLin : public FEInterpolation3d
     void surfaceEvaldNdxi(FloatMatrix &answer, int isurf, const FloatArray &lcoords) const override;
 
     void surfaceLocal2global(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const override;
-    void surfaceLocal2fullLocal( FloatArray &answer, int iSurf, const FloatArray &surfacelcoords, const FEICellGeometry &cellgeo ) const override;
     double surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const override;
     double surfaceGiveTransformationJacobian(int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const override;
     IntArray computeLocalSurfaceMapping(int iedge) const override;

src/oofemlib/feinterpol.h

@@ -429,18 +429,6 @@ class OOFEM_EXPORT FEInterpolation
      * @param cellgeo Underlying cell geometry.
      */
     virtual void boundaryLocal2Global(FloatArray &answer, int boundary, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const = 0;
-
-    /**
-    * Maps the local coordinates on boundary to local coordinates of entire geometry
-    * @param answer Full local coordinates.
-    * @param boundary Boundary number.
-    * @param boundarylcoords The local coordinates in the boundary local coordinate system
-    * @param cellgeo Underlying cell geometry.
-    */
-    virtual void boundaryLocal2fullLocal( FloatArray &answer, int boundary, const FloatArray &boundarylcoords, const FEICellGeometry &cellgeo ) const
-    {
-        OOFEM_ERROR( "Not implemented." );
-    }
 
     /**
      * Computes the integral @f$ \int_S n \cdot x \mathrm{d}s @f$.

src/oofemlib/feinterpol3d.C

@@ -90,11 +90,6 @@ void FEInterpolation3d :: boundaryLocal2Global(FloatArray &answer, int boundary,
     return this->surfaceLocal2global(answer, boundary, lcoords, cellgeo);
 }
 
-void FEInterpolation3d::boundaryLocal2fullLocal( FloatArray &answer, int boundary, const FloatArray &boundarylcoords, const FEICellGeometry &cellgeo ) const
-{
-    return this->surfaceLocal2fullLocal( answer, boundary, boundarylcoords, cellgeo );
-}
-
 IntArray FEInterpolation3d :: computeEdgeMapping(const IntArray &elemNodes, int iedge) const
 {
     const auto &ln = this->computeLocalEdgeMapping(iedge);

src/oofemlib/feinterpol3d.h

@@ -79,7 +79,6 @@ class OOFEM_EXPORT FEInterpolation3d : public FEInterpolation
     double boundaryEvalNormal(FloatArray &answer, int boundary, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const override;
     double boundaryGiveTransformationJacobian(int boundary, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const override;
     void boundaryLocal2Global(FloatArray &answer, int boundary, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const override;
-    void boundaryLocal2fullLocal( FloatArray &answer, int boundary, const FloatArray &boundarylcoords, const FEICellGeometry &cellgeo ) const override;
 
     /**@name Edge interpolation services */
     //@{
@@ -180,18 +179,6 @@ class OOFEM_EXPORT FEInterpolation3d : public FEInterpolation
      */
     virtual void surfaceLocal2global(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const = 0;
 
-    /**
-    * Maps the local coordinates on boundary to local coordinates of entire geometry
-    * @param answer Full local coordinates.
-    * @param boundary Boundary number.
-    * @param boundarylcoords The local coordinates in the boundary local coordinate system
-    * @param cellgeo Underlying cell geometry.
-    */
-    virtual void surfaceLocal2fullLocal( FloatArray &answer, int iSurf, const FloatArray &surfacelcoords, const FEICellGeometry &cellgeo ) const
-    {
-        OOFEM_ERROR( "Not implemented." );
-    }
-
     /**
      * Evaluates the edge jacobian of transformation between local and global coordinates.
      * @param isurf Determines the surface number.

src/sm/Elements/nlstructuralelement.C

@@ -64,10 +64,14 @@ void NLStructuralElement::computeDeformationGradientVectorAtBoundary( FloatArray
     // the receiver at time step tStep.
     // Order of components: 11, 22, 33, 23, 13, 12, 32, 31, 21 in the 3D.
 
+    IntArray bNodes = this->giveInterpolation()->boundaryGiveNodes( iSurf );
+
     // Obtain the current displacement vector of the element and subtract initial displacements (if present)
+
     FloatArray u;
-    this->computeVectorOf({ D_u, D_v, D_w }, VM_Total, tStep, u); // solution vector
+    this->computeBoundaryVectorOf(bNodes, { D_u, D_v, D_w }, VM_Total, tStep, u); // solution vector
     if ( initialDisplacements ) {
+        //possible source of future errors? Can this be subtracted from boundary u?
         u.subtract(* initialDisplacements);
     }
 

src/sm/Elements/structural3delement.C

@@ -80,14 +80,15 @@ Structural3DElement::computeBmatrixAt(GaussPoint *gp, FloatMatrix &answer, int l
 }
 
 void Structural3DElement::computeBHmatrixAtBoundary( GaussPoint *gp, FloatMatrix &answer, int iBoundary )
-// Does the same as computeBHmatrixAt, with the exception of having to calculate the missing coordinate of the boundary GP
+// Does the same as computeBHmatrixAt, with the exception of having to calculate the missing coordinate of the boundary GP and experiencing other nodes
 {
     FEInterpolation *interp = this->giveInterpolation();
     FloatMatrix dNdx;
-    FloatArray nCoords;
-    interp->boundaryLocal2fullLocal( nCoords, iBoundary, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
+    //FloatArray nCoords;
+    //interp->boundaryLocal2fullLocal( nCoords, iBoundary, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
 
-    interp->evaldNdx(dNdx, nCoords , FEIElementGeometryWrapper(this) );
+    interp->boundarySurfaceEvaldNdx( dNdx, iBoundary, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper( this ) );
+    //interp->evaldNdx( dNdx, nCoords, FEIElementGeometryWrapper( this ) );
 
     answer.resize(9, dNdx.giveNumberOfRows() * 3);
     answer.zero();

src/sm/Loads/hardmagneticboundarycondition.C

@@ -70,7 +70,9 @@ REGISTER_BoundaryCondition( HardMagneticBoundaryCondition );
         }
 
         FloatMatrix Ke;
-        IntArray r_loc, c_loc, bNodes;
+        IntArray r_loc, c_loc, bNodes, dofIDs;
+
+        dofIDs = domain->giveDefaultNodeDofIDArry();
 
         Set *set                   = this->giveDomain()->giveSet( this->set );
         const IntArray &boundaries = set->giveBoundaryList();
@@ -84,12 +86,27 @@ REGISTER_BoundaryCondition( HardMagneticBoundaryCondition );
                 bNodes = e->giveInterpolation()->boundarySurfaceGiveNodes( boundary );
             }
 
-            e->giveBoundaryLocationArray( r_loc, bNodes, this->dofs, r_s );
-            e->giveBoundaryLocationArray( c_loc, bNodes, this->dofs, c_s );
+            e->giveBoundaryLocationArray( r_loc, bNodes, dofIDs /* this->dofs*/, r_s );
+            e->giveBoundaryLocationArray( c_loc, bNodes, dofIDs /* this->dofs*/, c_s );
             this->computeTangentFromElement( Ke, e, boundary, tStep );
             answer.assemble( r_loc, c_loc, Ke );
+
+            //debug
+            if ( pos == 1 ) {
+                //compute numerical tangent
+                FloatMatrix Ke_num;
+                double pert = 1.e-6;
+                this->computeNumericalTangentFromElement( Ke_num, e, boundary, tStep, pert);
+
+                OOFEM_LOG_INFO( "Debugging surface %i\n", boundary );
+                OOFEM_LOG_INFO( "-------Analytical stiffness-----------\n" );
+                Ke.printYourself();
+                OOFEM_LOG_INFO( "-------Numerical stiffness------------\n" );
+                Ke_num.printYourself();
+
+            }
         }
-    }
+    }    
 
     void HardMagneticBoundaryCondition::assembleVector(FloatArray& answer , TimeStep* tStep ,
         CharType type , ValueModeType mode ,
@@ -278,5 +295,36 @@ REGISTER_BoundaryCondition( HardMagneticBoundaryCondition );
 
     }
 
-
-}//end namespace oofem;
+
+    //----------------------------------------------------------------------DEBUG-------------------------------------------------------------------------------------
+
+    void HardMagneticBoundaryCondition::computeNumericalTangentFromElement( FloatMatrix &answer, Element *e, int side, TimeStep *tStep, double perturb )
+    {
+        // debugging, numerically computing stiffness
+        answer.resize( 12, 12 );
+
+        FloatArray perturbedForceFront, perturbedForceBack;
+        for ( int i = 1; i <= 12; i++ ) {
+
+            this->computePerturbedLoadVectorFromElement( perturbedForceFront, e, side, tStep, perturb, i );
+            this->computePerturbedLoadVectorFromElement( perturbedForceBack, e, side, tStep, -perturb, i );
+
+            for ( int j = 1; j <= 12; j++ ) {
+
+                answer.at( i, j ) = ( perturbedForceFront.at( j ) - perturbedForceBack.at( j ) ) / 2 * perturb;
+            }
+        }
+    }
+
+    void HardMagneticBoundaryCondition::computePerturbedLoadVectorFromElement( FloatArray &answer, Element *e, int side, TimeStep *tStep, double perturb, int index )
+    {
+        // debugging, numerically perturbed load vector
+        IntArray boundaryNodes = e->giveInterpolation()->boundarySurfaceGiveNodes( side );
+
+        int ndof = 3 * boundaryNodes.giveSize();
+        answer.resize( ndof );
+
+        //todo
+    }
+
+    }//end namespace oofem;

src/sm/Loads/hardmagneticboundarycondition.h

@@ -140,5 +140,8 @@ class HardMagneticBoundaryCondition : public ActiveBoundaryCondition {
 private:
 
     void evaluateFreeSpaceStress();
+
+    void computeNumericalTangentFromElement( FloatMatrix &answer, Element *e, int side, TimeStep *tStep, double perturb ); //for debugging
+    void computePerturbedLoadVectorFromElement( FloatArray &answer, Element *e, int side, TimeStep *tStep, double perturb, int index);
 };
 } // end namespace oofem