diff --git a/src/oofemlib/fei3dhexalin.C b/src/oofemlib/fei3dhexalin.C index 651e294c5..af8171fe4 100644 --- a/src/oofemlib/fei3dhexalin.C +++ b/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 { diff --git a/src/oofemlib/fei3dhexalin.h b/src/oofemlib/fei3dhexalin.h index ad11f7c27..a2a60d9a1 100644 --- a/src/oofemlib/fei3dhexalin.h +++ b/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; diff --git a/src/oofemlib/feinterpol.h b/src/oofemlib/feinterpol.h index 919b27d98..a335843f4 100644 --- a/src/oofemlib/feinterpol.h +++ b/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$. diff --git a/src/oofemlib/feinterpol3d.C b/src/oofemlib/feinterpol3d.C index b2c97d575..eb2ded8a4 100644 --- a/src/oofemlib/feinterpol3d.C +++ b/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); diff --git a/src/oofemlib/feinterpol3d.h b/src/oofemlib/feinterpol3d.h index abae88e1d..e9c0a9cef 100644 --- a/src/oofemlib/feinterpol3d.h +++ b/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. diff --git a/src/sm/Elements/nlstructuralelement.C b/src/sm/Elements/nlstructuralelement.C index a22adf33f..0fca71633 100644 --- a/src/sm/Elements/nlstructuralelement.C +++ b/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); } diff --git a/src/sm/Elements/structural3delement.C b/src/sm/Elements/structural3delement.C index 2cdf104f4..fbf369348 100644 --- a/src/sm/Elements/structural3delement.C +++ b/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(); diff --git a/src/sm/Loads/hardmagneticboundarycondition.C b/src/sm/Loads/hardmagneticboundarycondition.C index 1f8e714df..7100a0f99 100644 --- a/src/sm/Loads/hardmagneticboundarycondition.C +++ b/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; diff --git a/src/sm/Loads/hardmagneticboundarycondition.h b/src/sm/Loads/hardmagneticboundarycondition.h index 79440a70a..1ce8bc9e9 100644 --- a/src/sm/Loads/hardmagneticboundarycondition.h +++ b/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