C++ FloatArray::plusProduct方法代码示例

void
SUPGElement2 :: computeDiffusionTerm_MB(FloatArray &answer, TimeStep *tStep)
{
    FloatArray u, eps, stress, bs, dDB_u;
    FloatMatrix b, un_gu, dDB;
    double Re = static_cast< FluidModel * >( domain->giveEngngModel() )->giveReynoldsNumber();

    answer.clear();

    this->computeVectorOfVelocities(VM_Total, tStep, u);

    int rule = 1;
    for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
        double dV = this->computeVolumeAround(gp);
        this->computeBMatrix(b, gp);
        this->computeDivTauMatrix(dDB, gp, tStep);
        this->computeUDotGradUMatrix( un_gu, gp, tStep->givePreviousStep() );
        eps.beProductOf(b, u);
        static_cast< FluidDynamicMaterial * >( this->giveMaterial() )->computeDeviatoricStressVector(stress, gp, eps, tStep);
        dDB_u.beProductOf(dDB, u);
        /* consistent part */
        answer.plusProduct(b, stress, dV / Re);

        /* SUPG term */
        answer.plusProduct(un_gu, dDB_u, ( -1.0 ) * t_supg * dV);
    }
}

void
SUPGElement2 :: computeBCRhsTerm_MB(FloatArray &answer, TimeStep *tStep)
{
    int nLoads;

    answer.clear();

    int rule = 0;
    IntegrationRule *iRule = this->integrationRulesArray [ rule ];
    FloatArray un, gVector, s, helpLoadVector;
    FloatMatrix b, nu;

    // add body load (gravity) termms
    nLoads = this->giveBodyLoadArray()->giveSize();
    for ( int i = 1; i <= nLoads; i++ ) {
        Load *load = domain->giveLoad( bodyLoadArray.at(i) );
        bcGeomType ltype = load->giveBCGeoType();
        if ( ( ltype == BodyLoadBGT ) && ( load->giveBCValType() == ForceLoadBVT ) ) {
            load->computeComponentArrayAt(gVector, tStep, VM_Total);
            if ( gVector.giveSize() ) {
                for ( GaussPoint *gp: *iRule ) {
                    this->computeUDotGradUMatrix( b, gp, tStep->givePreviousStep() );
                    this->computeNuMatrix(nu, gp);
                    double dV  = this->computeVolumeAround(gp);
                    double rho = this->giveMaterial()->give('d', gp);
                    answer.plusProduct(b, gVector, t_supg * rho * dV);
                    answer.plusProduct(nu, gVector, rho * dV);
                }
            }
        }
    }

    // integrate tractions
    // if no traction bc applied but side marked as with traction load
    // then zero traction is assumed !!!

    // loop over boundary load array
    nLoads = this->giveBoundaryLoadArray()->giveSize() / 2;
    for ( int i = 1; i <= nLoads; i++ ) {
        int n = boundaryLoadArray.at(1 + ( i - 1 ) * 2);
        int id = boundaryLoadArray.at(i * 2);
        Load *load  = domain->giveLoad(n);
        bcGeomType ltype = load->giveBCGeoType();
        if ( ltype == EdgeLoadBGT ) {
            this->computeEdgeLoadVector_MB(helpLoadVector, load, id, tStep);
            if ( helpLoadVector.giveSize() ) {
                answer.add(helpLoadVector);
            }
        } else if ( ltype == SurfaceLoadBGT ) {
            this->computeSurfaceLoadVector_MB(helpLoadVector, load, id, tStep);
            if ( helpLoadVector.giveSize() ) {
                answer.add(helpLoadVector);
            }
        } else {
            OOFEM_ERROR("unsupported load type class");
        }
    }
}

void MixedGradientPressureWeakPeriodic :: integrateTractionDev(FloatArray &answer, Element *el, int boundary, const FloatMatrix &ddev)
{
    // Computes the integral: int dt . dx dA
    FloatMatrix mMatrix;
    FloatArray normal, coords, vM_dev;

    FEInterpolation *interp = el->giveInterpolation(); // Geometry interpolation. The displacements or velocities must have the same interpolation scheme (on the boundary at least).

    int maxorder = this->order + interp->giveInterpolationOrder() * 3;
    std :: unique_ptr< IntegrationRule >ir( interp->giveBoundaryIntegrationRule(maxorder, boundary) );
    answer.clear();

    for ( GaussPoint *gp: *ir ) {
        const FloatArray &lcoords = gp->giveNaturalCoordinates();
        FEIElementGeometryWrapper cellgeo(el);

        double detJ = interp->boundaryEvalNormal(normal, boundary, lcoords, cellgeo);
        // Compute v_m = d_dev . x
        interp->boundaryLocal2Global(coords, boundary, lcoords, cellgeo);
        vM_dev.beProductOf(ddev, coords);

        this->constructMMatrix(mMatrix, coords, normal);

        answer.plusProduct(mMatrix, vM_dev, detJ * gp->giveWeight());
    }
}

void
StructuralInterfaceElementPhF :: giveInternalForcesVector_u(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
{

    
    FloatMatrix N;
    FloatArray u, traction, jump;
    
    IntArray dofIdArray;
    this->giveDofManDofIDMask_u(dofIdArray);
    this->computeVectorOf(dofIdArray, VM_Total, tStep, u);

    // subtract initial displacements, if defined
    if ( initialDisplacements.giveSize() ) {
        u.subtract(initialDisplacements);
    }
    
    // zero answer will resize accordingly when adding first contribution
    answer.clear();
    
    for ( auto &ip: *this->giveDefaultIntegrationRulePtr() ) {
        
        this->computeNmatrixAt(ip, N);
        
        jump.beProductOf(N, u);
        this->computeTraction(traction, ip, jump, tStep);
        //traction.resize(2);
        
        // compute internal cohesive forces as f = N^T*traction dA
        double dA = this->computeAreaAround(ip);
        answer.plusProduct(N, traction, dA);
    }

}

void
PhaseFieldElement :: giveInternalForcesVector_d(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
{
    // Computes int_V ( N^t *Nstress_d  +  B^t * g_c*l*grad(d)  ) dV
    FloatArray NStress, BStress, NS, a_d, grad_d;
    FloatMatrix N, B;
    NLStructuralElement *el = this->giveElement();
    this->computeDamageUnknowns( a_d, VM_Total, tStep );

    for ( auto &gp: el->giveIntegrationRule(0) ) {
        double dV = el->computeVolumeAround(gp);
        
        // compute generalized stress measures
        this->computeNd_matrixAt( *gp->giveNaturalCoordinates(), N);
        computeNStress_d(NStress, gp, tStep, useUpdatedGpRecord);
        NS.beTProductOf(N, NStress);
        answer.add(dV, NS);

        this->computeBd_matrixAt(gp, B);
        grad_d.beProductOf(B, a_d);
        double l = this->giveInternalLength();
        double g_c = this->giveCriticalEnergy();
        BStress = grad_d * l * g_c;
        answer.plusProduct(B, BStress, dV);
    }
}

void
MisesMatNl :: giveRemoteNonlocalStiffnessContribution(GaussPoint *gp, IntArray &rloc, const UnknownNumberingScheme &s,
                                                      FloatArray &rcontrib, TimeStep *tStep)
{
    double kappa, tempKappa;
    MisesMatNlStatus *status = static_cast< MisesMatNlStatus * >( this->giveStatus(gp) );
    StructuralElement *elem = static_cast< StructuralElement * >( gp->giveElement() );
    FloatMatrix b;
    LinearElasticMaterial *lmat = this->giveLinearElasticMaterial();
    double E = lmat->give('E', gp);

    elem->giveLocationArray(rloc, s);
    elem->computeBmatrixAt(gp, b);
    kappa = status->giveCumulativePlasticStrain();
    tempKappa = status->giveTempCumulativePlasticStrain();

    rcontrib.clear();
    if ( ( tempKappa - kappa ) > 0 ) {
        const FloatArray &stress = status->giveTempEffectiveStress();
        if ( gp->giveMaterialMode() == _1dMat ) {
            double coeff = sgn( stress.at(1) ) * E / ( E + H );
            rcontrib.plusProduct(b, stress, coeff);
            return;
        }
    }
    rcontrib.resize(b.giveNumberOfColumns());
    rcontrib.zero();
}

void
SUPGElement2 :: computeBCRhsTerm_MC(FloatArray &answer, TimeStep *tStep)
{
    int nLoads;
    FloatArray s, gVector, helpLoadVector;
    FloatMatrix g;

    int rule = 1;

    answer.clear();

    nLoads = this->giveBodyLoadArray()->giveSize();
    for ( int i = 1; i <= nLoads; i++ ) {
        Load *load  = domain->giveLoad( bodyLoadArray.at(i) );
        bcGeomType ltype = load->giveBCGeoType();
        if ( ( ltype == BodyLoadBGT ) && ( load->giveBCValType() == ForceLoadBVT ) ) {
            load->computeComponentArrayAt(gVector, tStep, VM_Total);
            if ( gVector.giveSize() ) {
                for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
                    this->computeGradPMatrix(g, gp);
                    double dV = this->computeVolumeAround(gp);
                    answer.plusProduct(g, gVector, t_pspg * dV);
                }
            }
        }
    }

    // integrate tractions
    // if no traction bc applied but side marked as with traction load
    // then zero traction is assumed !!!

    // loop over boundary load array
    nLoads = this->giveBoundaryLoadArray()->giveSize() / 2;
    for ( int i = 1; i <= nLoads; i++ ) {
        int n = boundaryLoadArray.at(1 + ( i - 1 ) * 2);
        int id = boundaryLoadArray.at(i * 2);
        Load *load = domain->giveLoad(n);
        bcGeomType ltype = load->giveBCGeoType();
        if ( ltype == EdgeLoadBGT ) {
            this->computeEdgeLoadVector_MC(helpLoadVector, load, id, tStep);
            if ( helpLoadVector.giveSize() ) {
                answer.add(helpLoadVector);
            }
        } else if ( ltype == SurfaceLoadBGT ) {
            this->computeSurfaceLoadVector_MC(helpLoadVector, load, id, tStep);
            if ( helpLoadVector.giveSize() ) {
                answer.add(helpLoadVector);
            }
        } else {
            OOFEM_ERROR("unsupported load type class");
        }
    }
}

void
Beam2d :: computeBoundaryEdgeLoadVector(FloatArray &answer, BoundaryLoad *load, int edge, CharType type, ValueModeType mode, TimeStep *tStep, bool global)
{
    answer.clear();

    if ( edge != 1 ) {
        OOFEM_ERROR("Beam2D only has 1 edge (the midline) that supports loads. Attempted to apply load to edge %d", edge);
    }

    if ( type != ExternalForcesVector ) {
        return;
    }

    double l = this->computeLength();
    FloatArray coords, t;
    FloatMatrix N, T;

    answer.clear();
    for ( auto &gp : *this->giveDefaultIntegrationRulePtr() ) {
        const FloatArray &lcoords = gp->giveNaturalCoordinates();
        this->computeNmatrixAt(lcoords, N);
        if ( load ) {
            this->computeGlobalCoordinates(coords, lcoords);
            load->computeValues(t, tStep, coords, { D_u, D_w, R_v }, mode);
        } else {
            load->computeValues(t, tStep, lcoords, { D_u, D_w, R_v }, mode);
        }

        if ( load->giveCoordSystMode() == Load :: CST_Global ) {
            if ( this->computeLoadGToLRotationMtrx(T) ) {
                t.rotatedWith(T, 'n');
            }
        }

        double dl = gp->giveWeight() * 0.5 * l;
        answer.plusProduct(N, t, dl);
    }

    if (global) {
      // Loads from sets expects global c.s.
      this->computeGtoLRotationMatrix(T);
      answer.rotatedWith(T, 't');
    }
}

void
StructuralInterfaceElement :: giveInternalForcesVector(FloatArray &answer,
                                                       TimeStep *tStep, int useUpdatedGpRecord)
{
    // Computes internal forces
    // if useGpRecord == 1 then data stored in ip->giveStressVector() are used
    // instead computing stressVector through this->ComputeStressVector();
    // this must be done after you want internal forces after element->updateYourself()
    // has been called for the same time step.

    IntegrationRule *iRule = integrationRulesArray [ giveDefaultIntegrationRule() ];

    FloatMatrix N, rotationMatGtoL;
    FloatArray u, traction, tractionTemp, jump;

    this->computeVectorOf(VM_Total, tStep, u);
    // subtract initial displacements, if defined
    if ( initialDisplacements ) {
        u.subtract(* initialDisplacements);
    }

    // zero answer will resize accordingly when adding first contribution
    answer.clear();

    for ( IntegrationPoint *ip: *iRule ) {
        this->computeNmatrixAt(ip, N);
        //if ( useUpdatedGpRecord == 1 ) {
        //    StructuralInterfaceMaterialStatus *status = static_cast< StructuralInterfaceMaterialStatus * >( ip->giveMaterialStatus() );
        //    //temp
        //    FloatArray traction3d;
        //    traction3d = status->giveTraction();
        //    traction.setValues(2, traction3d.at(1), traction3d.at(3) );
        //} else {
            jump.beProductOf(N, u);
            this->computeTraction(traction, ip, jump, tStep);
        //}

        // compute internal cohesive forces as f = N^T*traction dA
        double dA = this->computeAreaAround(ip);
        answer.plusProduct(N, traction, dA);
    }

}

int
MisesMatNl :: giveLocalNonlocalStiffnessContribution(GaussPoint *gp, IntArray &loc, const UnknownNumberingScheme &s,
                                                     FloatArray &lcontrib, TimeStep *tStep)
{
    double nlKappa, damage, tempDamage, dDamF;
    MisesMatNlStatus *status = static_cast< MisesMatNlStatus * >( this->giveStatus(gp) );
    StructuralElement *elem = static_cast< StructuralElement * >( gp->giveElement() );
    FloatMatrix b;

    this->computeCumPlasticStrain(nlKappa, gp, tStep);
    damage = status->giveDamage();
    tempDamage = status->giveTempDamage();
    if ( ( tempDamage - damage ) > 0 ) {
        const FloatArray &stress = status->giveTempEffectiveStress();
        elem->giveLocationArray(loc, s);
        elem->computeBmatrixAt(gp, b);
        dDamF = computeDamageParamPrime(nlKappa);
        lcontrib.clear();
        lcontrib.plusProduct(b, stress, mm * dDamF);
    }

    return 1;
}

void
NLStructuralElement :: giveInternalForcesVector_withIRulesAsSubcells(FloatArray &answer,
                                                                     TimeStep *tStep, int useUpdatedGpRecord)
{
    /*
     * Returns nodal representation of real internal forces computed from first Piola-Kirchoff stress
     * if useGpRecord == 1 then stresses stored in the gp are used, otherwise stresses are computed
     * this must be done if you want internal forces after element->updateYourself() has been called
     * for the same time step.
     * The integration procedure uses an integrationRulesArray for numerical integration.
     * Each integration rule is considered to represent a separate sub-cell/element. Typically this would be used when
     * integration of the element domain needs special treatment, e.g. when using the XFEM.
     */

    FloatMatrix B;
    FloatArray vStress, vStrain;

    IntArray irlocnum;
    FloatArray *m = & answer, temp;
    if ( this->giveInterpolation() && this->giveInterpolation()->hasSubPatchFormulation() ) {
        m = & temp;
    }

    // zero answer will resize accordingly when adding first contribution
    answer.clear();


    // loop over individual integration rules
    for ( auto &iRule: integrationRulesArray ) {
        for ( GaussPoint *gp: *iRule ) {
            StructuralMaterialStatus *matStat = static_cast< StructuralMaterialStatus * >( gp->giveMaterialStatus() );

            if ( nlGeometry == 0 ) {
                this->computeBmatrixAt(gp, B);
                if ( useUpdatedGpRecord == 1 ) {
                    vStress = matStat->giveStressVector();
                } else {
                    this->computeStrainVector(vStrain, gp, tStep);
                    this->computeStressVector(vStress, vStrain, gp, tStep);
                }
            } else if ( nlGeometry == 1 ) {
                if ( this->domain->giveEngngModel()->giveFormulation() == AL ) { // Cauchy stress
                    if ( useUpdatedGpRecord == 1 ) {
                        vStress = matStat->giveCVector();
                    } else {
                        this->computeCauchyStressVector(vStress, gp, tStep);
                    }

                    this->computeBmatrixAt(gp, B);
                } else { // First Piola-Kirchhoff stress
                    if ( useUpdatedGpRecord == 1 ) {
                        vStress = matStat->givePVector();
                    } else {
                        this->computeFirstPKStressVector(vStress, gp, tStep);
                    }

                    this->computeBHmatrixAt(gp, B);
                }
            }

            if ( vStress.giveSize() == 0 ) { //@todo is this really necessary?
                break;
            }

            // compute nodal representation of internal forces at nodes as f = B^T*stress dV
            double dV = this->computeVolumeAround(gp);
            m->plusProduct(B, vStress, dV);

            // localize irule contribution into element matrix
            if ( this->giveIntegrationRuleLocalCodeNumbers(irlocnum, *iRule) ) {
                answer.assemble(* m, irlocnum);
                m->clear();
            }
        }
    }

    // if inactive: update fields but do not give any contribution to the structure
    if ( !this->isActivated(tStep) ) {
        answer.zero();
        return;
    }
}

void
NLStructuralElement :: giveInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
{
    FloatMatrix B;
    FloatArray vStress, vStrain, u;

    // This function can be quite costly to do inside the loops when one has many slave dofs.
    this->computeVectorOf(VM_Total, tStep, u);
    // subtract initial displacements, if defined
    if ( initialDisplacements ) {
        u.subtract(* initialDisplacements);
    }

    // zero answer will resize accordingly when adding first contribution
    answer.clear();

    for ( auto &gp: *this->giveDefaultIntegrationRulePtr() ) {
        StructuralMaterialStatus *matStat = static_cast< StructuralMaterialStatus * >( gp->giveMaterialStatus() );

        // Engineering (small strain) stress
        if ( nlGeometry == 0 ) {
            this->computeBmatrixAt(gp, B);
            if ( useUpdatedGpRecord == 1 ) {
                vStress = matStat->giveStressVector();
            } else {
                ///@todo Is this really what we should do for inactive elements?
                if ( !this->isActivated(tStep) ) {
                    vStrain.resize( StructuralMaterial :: giveSizeOfVoigtSymVector( gp->giveMaterialMode() ) );
                    vStrain.zero();
                }
                vStrain.beProductOf(B, u);
                this->computeStressVector(vStress, vStrain, gp, tStep);
            }
        } else if ( nlGeometry == 1 ) {  // First Piola-Kirchhoff stress
            if ( this->domain->giveEngngModel()->giveFormulation() == AL ) { // Cauchy stress
                if ( useUpdatedGpRecord == 1 ) {
                    vStress = matStat->giveCVector();
                } else {
                    this->computeCauchyStressVector(vStress, gp, tStep);
                }

                this->computeBmatrixAt(gp, B);
            } else { // First Piola-Kirchhoff stress
                if ( useUpdatedGpRecord == 1 ) {
                    vStress = matStat->givePVector();
                } else {
                    this->computeFirstPKStressVector(vStress, gp, tStep);
                    ///@todo This is actaully inefficient since it constructs B and twice and collects the nodal unknowns over and over.
                }

                this->computeBHmatrixAt(gp, B);
            }
        }

        if ( vStress.giveSize() == 0 ) { /// @todo is this check really necessary?
            break;
        }

        // Compute nodal internal forces at nodes as f = B^T*Stress dV
        double dV  = this->computeVolumeAround(gp);

        if ( nlGeometry == 1 ) {  // First Piola-Kirchhoff stress
            if ( vStress.giveSize() == 9 ) {
                FloatArray stressTemp;
                StructuralMaterial :: giveReducedVectorForm( stressTemp, vStress, gp->giveMaterialMode() );
                answer.plusProduct(B, stressTemp, dV);
            } else   {
                answer.plusProduct(B, vStress, dV);
            }
        } else   {
            if ( vStress.giveSize() == 6 ) {
                // It may happen that e.g. plane strain is computed
                // using the default 3D implementation. If so,
                // the stress needs to be reduced.
                // (Note that no reduction will take place if
                //  the simulation is actually 3D.)
                FloatArray stressTemp;
                StructuralMaterial :: giveReducedSymVectorForm( stressTemp, vStress, gp->giveMaterialMode() );
                answer.plusProduct(B, stressTemp, dV);
            } else   {
                answer.plusProduct(B, vStress, dV);
            }
        }
    }

    // If inactive: update fields but do not give any contribution to the internal forces
    if ( !this->isActivated(tStep) ) {
        answer.zero();
        return;
    }
}

//.........这里部分代码省略.........
        OOFEM_ERROR("No interpolation or default integration available for element.");
    }
    std :: unique_ptr< IntegrationRule > iRule( fei->giveBoundaryIntegrationRule(fei->giveInterpolationOrder()-1, side) );

    int nsd = e->giveDomain()->giveNumberOfSpatialDimensions();

    if ( side == -1 ) {
        side = 1;
    }

    answer.clear();

    if ( nsd == 2 ) {

        FEInterpolation2d *fei2d = static_cast< FEInterpolation2d * >(fei);

        ///@todo More of this grunt work should be moved to the interpolation classes
        IntArray bnodes;
        fei2d->boundaryGiveNodes(bnodes, side);
        int nodes = bnodes.giveSize();
        FloatMatrix xy(2, nodes);
        for ( int i = 1; i <= nodes; i++ ) {
            Node *node = e->giveNode(bnodes.at(i));
            ///@todo This should rely on the xindex and yindex in the interpolator..
            xy.at(1, i) = node->giveCoordinate(1);
            xy.at(2, i) = node->giveCoordinate(2);
        }

        FloatArray tmp; // Integrand
        FloatArray es; // Tangent vector to curve
        FloatArray dNds;

        if ( e->giveDomain()->isAxisymmetric() ) {
            FloatArray N;
            FloatArray gcoords;
            for ( GaussPoint *gp: *iRule ) {
                fei2d->edgeEvaldNds( dNds, side, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e) );
                fei->boundaryEvalN( N, side, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e) );
                double J = fei->boundaryGiveTransformationJacobian( side, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e) );
                fei->boundaryLocal2Global( gcoords, side, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e) );
                double r = gcoords(0); // First coordinate is the radial coord.

                es.beProductOf(xy, dNds);

                tmp.resize( 2 * nodes);
                for ( int i = 0; i < nodes; i++ ) {
                    tmp(2 * i)   = dNds(i) * es(0) * r + N(i);
                    tmp(2 * i + 1) = dNds(i) * es(1) * r;
                }

                answer.add(- 2 * M_PI * gamma * J * gp->giveWeight(), tmp);
            }
        } else {
            for ( GaussPoint *gp: *iRule ) {
                double t = e->giveCrossSection()->give(CS_Thickness, gp);
                fei2d->edgeEvaldNds( dNds, side, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e) );
                double J = fei->boundaryGiveTransformationJacobian( side, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e) );
                es.beProductOf(xy, dNds);

                tmp.resize( 2 * nodes);
                for ( int i = 0; i < nodes; i++ ) {
                    tmp(2 * i)   = dNds(i) * es(0);
                    tmp(2 * i + 1) = dNds(i) * es(1);
                    //B.at(1, 1+i*2) = B.at(2, 2+i*2) = dNds(i);
                }
                //tmp.beTProductOf(B, es);

                answer.add(- t * gamma * J * gp->giveWeight(), tmp);
            }
        }
    } else if ( nsd ==  3 ) {

        FEInterpolation3d *fei3d = static_cast< FEInterpolation3d * >(fei);
        FloatArray n, surfProj;
        FloatMatrix dNdx, B;
        for ( GaussPoint *gp: *iRule ) {
            fei3d->surfaceEvaldNdx( dNdx, side, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e) );
            double J = fei->boundaryEvalNormal( n, side, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e) );

            // [I - n(x)n]  in voigt form:
            surfProj = {1. - n(0)*n(0), 1. - n(1)*n(1), 1. - n(2)*n(2),
                        - n(1)*n(2), - n(0)*n(2), - n(0)*n(1),
            };

            // Construct B matrix of the surface nodes
            B.resize(6, 3 * dNdx.giveNumberOfRows());
            for ( int i = 1; i <= dNdx.giveNumberOfRows(); i++ ) {
                B.at(1, 3 * i - 2) = dNdx.at(i, 1);
                B.at(2, 3 * i - 1) = dNdx.at(i, 2);
                B.at(3, 3 * i - 0) = dNdx.at(i, 3);

                B.at(5, 3 * i - 2) = B.at(4, 3 * i - 1) = dNdx.at(i, 3);
                B.at(6, 3 * i - 2) = B.at(4, 3 * i - 0) = dNdx.at(i, 2);
                B.at(6, 3 * i - 1) = B.at(5, 3 * i - 0) = dNdx.at(i, 1);
            }

            answer.plusProduct(B, surfProj, -gamma * J * gp->giveWeight() );
        }
    }
}

void
IntElLine2IntPen :: giveInternalForcesVector(FloatArray &answer,
                                                       TimeStep *tStep, int useUpdatedGpRecord)
{
    // Computes internal forces
	// For this element we use an "interior penalty" formulation, where
	// the cohesive zone contribution is weakened, i.e. the traction and
	// test function for the cohesive zone are projected onto a reduced
	// space. The purpose of the projection is to improve the stability
	// properties of the formulation, thereby avoiding traction oscilations.

    FloatMatrix N;
    FloatArray u, traction, jump;

    this->computeVectorOf(VM_Total, tStep, u);
    // subtract initial displacements, if defined
    if ( initialDisplacements.giveSize() ) {
        u.subtract(initialDisplacements);
    }

    // zero answer will resize accordingly when adding first contribution
    answer.clear();



    // First loop over GP: compute projection of test function and traction.
    // The setting is as follows: we have an interface with quadratic interpolation and we
    // wish to project onto the space of constant functions on each element.

    // The projection of t becomes a constant
//    FloatArray proj_t;
//    proj_t.clear();


    FloatArray proj_jump;
    proj_jump.clear();

    // Projecting the basis functions gives a constant for each basis function.
    FloatMatrix proj_N;
    proj_N.clear();

    double area = 0.;

    for ( auto &ip: *this->giveDefaultIntegrationRulePtr() ) {

        this->computeNmatrixAt(ip, N);
        jump.beProductOf(N, u);
//        this->computeTraction(traction, ip, jump, tStep);

        double dA = this->computeAreaAround(ip);
        area += dA;

        proj_jump.add(dA, jump);
        proj_N.add(dA, N);
    }

//    printf("area: %e\n", area);
    proj_jump.times(1./area);
    proj_N.times(1./area);

//    printf("proj_N: "); proj_N.printYourself();


    // Second loop over GP: assemble contribution to internal forces as we are used to.
    for ( auto &ip: *this->giveDefaultIntegrationRulePtr() ) {
//        this->computeNmatrixAt(ip, N);
//        jump.beProductOf(N, u);
        this->computeTraction(traction, ip, proj_jump, tStep);

        // compute internal cohesive forces as f = N^T*traction dA
        double dA = this->computeAreaAround(ip);
        answer.plusProduct(proj_N, traction, dA);
    }

}