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

void PrescribedGradientBCPeriodic :: computeField(FloatArray &sigma, TimeStep *tStep)
{
    DofIDEquationNumbering pnum(true, strain_id);
    EngngModel *emodel = this->giveDomain()->giveEngngModel();
    FloatArray tmp, sig_tmp;
    int npeq = strain_id.giveSize();
    // sigma = residual (since we use the slave dofs) = f_ext - f_int
    sig_tmp.resize(npeq);
    sig_tmp.zero();
    emodel->assembleVector(sig_tmp, tStep, InternalForceAssembler(), VM_Total, pnum, this->domain);
    tmp.resize(npeq);
    tmp.zero();
    emodel->assembleVector(tmp, tStep, ExternalForceAssembler(), VM_Total, pnum, this->domain);
    sig_tmp.subtract(tmp);
    // Divide by the RVE-volume
    sig_tmp.times(1.0 / ( this->domainSize(this->giveDomain(), this->set) + this->domainSize(this->giveDomain(), this->masterSet) ));

    sigma.resize(sig_tmp.giveSize());
    if ( sig_tmp.giveSize() == 9 ) {
        sigma.assemble(sig_tmp, {1, 9, 8, 6, 2, 7, 5, 4, 3});
    } else if ( sig_tmp.giveSize() == 4 ) {
        sigma.assemble(sig_tmp, {1, 4, 3, 2});
    } else {
        sigma = sig_tmp;
    }
}

void
Quad1MindlinShell3D :: giveInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
{
    // We need to overload this for practical reasons (this 3d shell has all 9 dofs, but the shell part only cares for the first 8)
    // This elements adds an additional stiffness for the so called drilling dofs, meaning we need to work with all 9 components.
    FloatMatrix b, d;
    FloatArray n, strain, stress;
    FloatArray shellUnknowns(20), drillUnknowns(4), unknowns;

    this->computeVectorOf(EID_MomentumBalance, VM_Total, tStep, unknowns);
    // Split this for practical reasons into normal shell dofs and drilling dofs
    for ( int i = 0; i < 4; ++i ) {
        shellUnknowns(0 + i*5) = unknowns(0 + i*6);
        shellUnknowns(1 + i*5) = unknowns(1 + i*6);
        shellUnknowns(2 + i*5) = unknowns(2 + i*6);
        shellUnknowns(3 + i*5) = unknowns(3 + i*6);
        shellUnknowns(4 + i*5) = unknowns(4 + i*6);
        drillUnknowns(i) = unknowns(5 + i*6);
    }

    FloatArray shellForces(20), drillMoment(4);
    shellForces.zero();
    drillMoment.zero();
    StructuralCrossSection *cs = this->giveStructuralCrossSection();
    double drillCoeff = cs->give(CS_DrillingStiffness);

    IntegrationRule *iRule = integrationRulesArray [ 0 ];
    for ( int i = 0; i < iRule->giveNumberOfIntegrationPoints(); i++ ) {
        GaussPoint *gp = iRule->getIntegrationPoint(i);
        this->computeBmatrixAt(gp, b);
        double dV = this->computeVolumeAround(gp);

        if ( useUpdatedGpRecord ) {
            stress = static_cast< StructuralMaterialStatus * >( this->giveMaterial()->giveStatus(gp) )->giveStressVector();
        } else {
            strain.beProductOf(b, shellUnknowns);
            cs->giveRealStress_Shell(stress, gp, strain, tStep);
        }
        shellForces.plusProduct(b, stress, dV);

        // Drilling stiffness is here for improved numerical properties
        if (drillCoeff > 0.) {
            this->interp.evalN(n, *gp->giveCoordinates(), FEIVoidCellGeometry());
            for ( int j = 0; j < 4; j++) {
                n(j) -= 0.25;
            }
            double dtheta = n.dotProduct(drillUnknowns);
            drillMoment.add(drillCoeff * dV * dtheta, n); ///@todo Decide on how to alpha should be defined.
        }
    }

    answer.resize(24);
    answer.zero();
    answer.assemble(shellForces, this->shellOrdering);

    if (drillCoeff > 0.) {
        answer.assemble(drillMoment, this->drillOrdering);
    }
}

void Tr21Stokes :: computeBodyLoadVectorAt(FloatArray &answer, Load *load, TimeStep *tStep)
{
    IntegrationRule *iRule = this->integrationRulesArray [ 0 ];
    GaussPoint *gp;
    FloatArray N, gVector, *lcoords, temparray(15);
    double dA, detJ, rho;

    load->computeComponentArrayAt(gVector, tStep, VM_Total);
    temparray.zero();
    if ( gVector.giveSize() ) {
        for ( int k = 0; k < iRule->getNumberOfIntegrationPoints(); k++ ) {
            gp = iRule->getIntegrationPoint(k);
            lcoords = gp->giveCoordinates();

            rho = this->giveMaterial()->giveCharacteristicValue(MRM_Density, gp, tStep);
            detJ = fabs( this->interpolation_quad.giveTransformationJacobian(* lcoords, FEIElementGeometryWrapper(this)) );
            dA = detJ * gp->giveWeight();

            this->interpolation_quad.evalN(N, * lcoords, FEIElementGeometryWrapper(this));
            for ( int j = 0; j < 6; j++ ) {
                temparray(2 * j)     += N(j) * rho * gVector(0) * dA;
                temparray(2 * j + 1) += N(j) * rho * gVector(1) * dA;
            }
        }
    }

    answer.resize(15);
    answer.zero();
    answer.assemble( temparray, this->ordering );
}

void
ContactDefinition :: computeContactForces(FloatArray &answer, TimeStep *tStep, CharType type, ValueModeType mode,
                                const UnknownNumberingScheme &s, Domain *domain, FloatArray *eNorms)
{
    //Loop through all the contact elements and let them return their internal forces vector
    FloatArray Fc;
    IntArray locArray;
    
    // TODO ask masters that are potentially in contact and not everyone
    for ( ContactElement *master : this->masterElementList ) {
        
        // These acts as external forces so move them to the lhs
        master->computeContactForces(Fc, tStep, type, mode, s, domain, eNorms);
        Fc.negated();
        
        if ( Fc.giveSize() ) {
            master->giveLocationArray(locArray, s);
            answer.assemble(Fc, locArray);
        
          if ( eNorms ) {
              eNorms->assembleSquared( Fc, locArray );
          }
          
        }
    }
  
  
}

void
NonStationaryTransportProblem :: assembleDirichletBcRhsVector(FloatArray &answer, TimeStep *tStep,
                                                              ValueModeType mode,
                                                              const UnknownNumberingScheme &ns, Domain *d)
{
    IntArray loc, dofids;
    FloatArray rp, charVec;
    FloatMatrix s;
    FloatMatrix capacity;

    int nelem = d->giveNumberOfElements();

    for ( int ielem = 1; ielem <= nelem; ielem++ ) {
        Element *element = d->giveElement(ielem);

        element->giveElementDofIDMask(dofids);
        element->computeVectorOfPrescribed(dofids, mode, tStep, rp);
        if ( rp.containsOnlyZeroes() ) {
            continue;
        } else {
            element->giveCharacteristicMatrix(s, TangentStiffnessMatrix, tStep);
            element->giveCharacteristicMatrix(capacity, lumpedCapacityStab ? LumpedMassMatrix : MassMatrix, tStep);
            s.times(this->alpha);
            s.add(1. / tStep->giveTimeIncrement(), capacity);

            charVec.beProductOf(s, rp);
            charVec.negated();

            element->giveLocationArray(loc, ns);
            answer.assemble(charVec, loc);
        }
    } // end element loop
}

void SurfaceTensionBoundaryCondition :: assembleVector(FloatArray &answer, TimeStep *tStep,
                                                       CharType type, ValueModeType mode,
                                                       const UnknownNumberingScheme &s, FloatArray *eNorms)
{
    if ( type != ExternalForcesVector ) {
        return;
    }

    FloatArray fe;
    IntArray loc, masterdofids, bNodes;

    Set *set = this->giveDomain()->giveSet(this->set);
    const IntArray &boundaries = set->giveBoundaryList();

    for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
        Element *e = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
        int boundary = boundaries.at(pos * 2);

        e->giveInterpolation()->boundaryGiveNodes(bNodes, boundary);

        e->giveBoundaryLocationArray(loc, bNodes, this->dofs, s, & masterdofids);
        this->computeLoadVectorFromElement(fe, e, boundary, tStep);
        answer.assemble(fe, loc);
        if ( eNorms ) {
            eNorms->assembleSquared(fe, masterdofids);
        }
    }
}

void
NonStationaryTransportProblem :: assembleDirichletBcRhsVector(FloatArray &answer, TimeStep *tStep, EquationID ut,
                                                              ValueModeType mode, CharType lhsType,
                                                              const UnknownNumberingScheme &ns, Domain *d)
{
    int ielem;
    IntArray loc;
    Element *element;
    FloatArray rp, charVec;
    FloatMatrix s, bcMtrx;

    int nelem = d->giveNumberOfElements();

    for ( ielem = 1; ielem <= nelem; ielem++ ) {
        element = d->giveElement(ielem);

        element->computeVectorOfPrescribed(EID_ConservationEquation, mode, tStep, rp);
        if ( rp.containsOnlyZeroes() ) {
            continue;
        } else {
            this->giveElementCharacteristicMatrix(s, ielem, lhsType, tStep, d);
            element->giveCharacteristicMatrix(bcMtrx, LHSBCMatrix, tStep);
            s.add(bcMtrx);
            charVec.beProductOf(s, rp);
            charVec.negated();

            element->giveLocationArray(loc, ut, ns);
            answer.assemble(charVec, loc);
        }
    } // end element loop

}

void Tr21Stokes :: computeLoadVector(FloatArray &answer, BodyLoad *load, CharType type, ValueModeType mode, TimeStep *tStep)
{
    if ( type != ExternalForcesVector ) {
        answer.clear();
        return;
    }

    FluidDynamicMaterial *mat = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
    FloatArray N, gVector, temparray(12);

    load->computeComponentArrayAt(gVector, tStep, VM_Total);
    temparray.zero();
    if ( gVector.giveSize() ) {
        for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
            const FloatArray &lcoords = gp->giveNaturalCoordinates();

            double rho = mat->give('d', gp);
            double detJ = fabs( this->interpolation_quad.giveTransformationJacobian( lcoords, FEIElementGeometryWrapper(this) ) );
            double dA = detJ * gp->giveWeight();

            this->interpolation_quad.evalN( N, lcoords, FEIElementGeometryWrapper(this) );
            for ( int j = 0; j < 6; j++ ) {
                temparray(2 * j)     += N(j) * rho * gVector(0) * dA;
                temparray(2 * j + 1) += N(j) * rho * gVector(1) * dA;
            }
        }
    }

    answer.resize(15);
    answer.zero();
    answer.assemble(temparray, this->momentum_ordering);
}

void
GradDpElement :: computeForceLoadVector(FloatArray &answer, TimeStep *tStep, ValueModeType mode)
{
    //set displacement and nonlocal location array
    this->setDisplacementLocationArray();
    this->setNonlocalLocationArray();


    FloatArray localForces(locSize);
    FloatArray nlForces(nlSize);
    answer.resize(totalSize);

    this->computeLocForceLoadVector(localForces, tStep, mode);

    answer.assemble(localForces, locU);
    answer.assemble(nlForces, locK);
}

void
TR_SHELL02 :: giveCharacteristicVector(FloatArray &answer, CharType mtrx, ValueModeType mode, TimeStep *tStep)
//
// returns characteristics vector of receiver accordind to mtrx
//
{
    FloatArray aux;

    answer.resize(18);
    answer.zero();

    plate->giveCharacteristicVector(aux, mtrx, mode, tStep);
    if ( !aux.isEmpty() ) answer.assemble(aux, loc_plate);

    membrane->giveCharacteristicVector(aux, mtrx, mode, tStep);
    if ( !aux.isEmpty() ) answer.assemble(aux, loc_membrane);
}

void
TR_SHELL01 :: giveCharacteristicVector(FloatArray &answer, CharType mtrx, ValueModeType mode, TimeStep *tStep)
//
// returns characteristics vector of receiver accordind to mtrx
//
{
  IntArray loc(9);
  FloatArray aux;

  answer.resize(18);
  answer.zero();

  plate->giveCharacteristicVector(aux, mtrx, mode, tStep);
  loc.setValues(9, 3,4,5, 9,10,11, 15,16,17);
  answer.assemble(aux, loc);

  membrane->giveCharacteristicVector(aux, mtrx, mode, tStep);
  loc.setValues(9, 1,2,6, 7,8,12, 13,14,18);
  answer.assemble(aux, loc);
}

void Tr21Stokes :: computeInternalForcesVector(FloatArray &answer, TimeStep *tStep)
{
    FluidDynamicMaterial *mat = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
    FloatArray a_pressure, a_velocity, devStress, epsp, Nh, dNv(12);
    double r_vol, pressure;
    FloatMatrix dN, B(3, 12);
    B.zero();

    this->computeVectorOfVelocities(VM_Total, tStep, a_velocity);
    this->computeVectorOfPressures(VM_Total, tStep, a_pressure);

    FloatArray momentum, conservation;

    for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
        const FloatArray &lcoords = gp->giveNaturalCoordinates();

        double detJ = fabs( this->interpolation_quad.evaldNdx( dN, lcoords, FEIElementGeometryWrapper(this) ) );
        this->interpolation_lin.evalN( Nh, lcoords, FEIElementGeometryWrapper(this) );
        double dA = detJ * gp->giveWeight();

        for ( int j = 0, k = 0; j < dN.giveNumberOfRows(); j++, k += 2 ) {
            dNv(k)     = B(0, k)     = B(2, k + 1) = dN(j, 0);
            dNv(k + 1) = B(1, k + 1) = B(2, k)     = dN(j, 1);
        }

        pressure = Nh.dotProduct(a_pressure);
        epsp.beProductOf(B, a_velocity);

        mat->computeDeviatoricStressVector(devStress, r_vol, gp, epsp, pressure, tStep);

        momentum.plusProduct(B, devStress, dA);
        momentum.add(-pressure * dA, dNv);
        conservation.add(r_vol * dA, Nh);
    }

    answer.resize(15);
    answer.zero();
    answer.assemble(momentum, this->momentum_ordering);
    answer.assemble(conservation, this->conservation_ordering);
}

void Tr21Stokes :: computeInternalForcesVector(FloatArray &answer, TimeStep *tStep)
{
    IntegrationRule *iRule = integrationRulesArray [ 0 ];
    FluidDynamicMaterial *mat = ( FluidDynamicMaterial * ) this->domain->giveMaterial(this->material);
    FloatArray a_pressure, a_velocity, devStress, epsp, BTs, Nh, dNv(12);
    double r_vol, pressure;
    FloatMatrix dN, B(3, 12);
    B.zero();
    
    this->computeVectorOf(EID_MomentumBalance, VM_Total, tStep, a_velocity);
    this->computeVectorOf(EID_ConservationEquation, VM_Total, tStep, a_pressure);
    
    FloatArray momentum(12), conservation(3);
    momentum.zero();
    conservation.zero();
    GaussPoint *gp;

    for ( int i = 0; i < iRule->getNumberOfIntegrationPoints(); i++ ) {
        gp = iRule->getIntegrationPoint(i);
        FloatArray *lcoords = gp->giveCoordinates();

        double detJ = fabs(this->interpolation_quad.giveTransformationJacobian(* lcoords, FEIElementGeometryWrapper(this)));
        this->interpolation_quad.evaldNdx(dN, * lcoords, FEIElementGeometryWrapper(this));
        this->interpolation_lin.evalN(Nh, * lcoords, FEIElementGeometryWrapper(this));
        double dA = detJ * gp->giveWeight();

        for ( int j = 0, k = 0; j < 6; j++, k += 2 ) {
            dNv(k)     = B(0, k)     = B(2, k + 1) = dN(j, 0);
            dNv(k + 1) = B(1, k + 1) = B(2, k)     = dN(j, 1);
        }

        pressure = Nh.dotProduct(a_pressure);
        epsp.beProductOf(B, a_velocity);

        mat->computeDeviatoricStressVector(devStress, r_vol, gp, epsp, pressure, tStep);
        BTs.beTProductOf(B, devStress);

        momentum.add(dA, BTs);
        momentum.add(-pressure*dA, dNv);
        conservation.add(r_vol*dA, Nh);
    }

    FloatArray temp(15);
    temp.zero();
    temp.addSubVector(momentum, 1);
    temp.addSubVector(conservation, 13);

    answer.resize(15);
    answer.zero();
    answer.assemble(temp, this->ordering);
}

void Tr1Darcy :: computeEdgeBCSubVectorAt(FloatArray &answer, Load *load, int iEdge, TimeStep *tStep)
{
    /*
     * Given the load *load, return it's contribution.
     *
     */

    answer.resize(3);
    answer.zero();

    if ( load->giveType() == TransmissionBC ) {                 // Neumann boundary conditions (traction)
        BoundaryLoad *boundaryLoad;
        boundaryLoad = ( BoundaryLoad * ) load;

        int numberOfEdgeIPs;
        numberOfEdgeIPs = ( int ) ceil( ( boundaryLoad->giveApproxOrder() + 1. ) / 2. ) * 2;

        GaussIntegrationRule iRule(1, this, 1, 1);
        GaussPoint *gp;
        FloatArray N, loadValue, reducedAnswer;
        reducedAnswer.resize(3);
        reducedAnswer.zero();
        IntArray mask;

        iRule.setUpIntegrationPoints(_Line, numberOfEdgeIPs, _Unknown);

        for ( int i = 0; i < iRule.getNumberOfIntegrationPoints(); i++ ) {
            gp = iRule.getIntegrationPoint(i);
            FloatArray *lcoords = gp->giveCoordinates();
            this->interpolation_lin.edgeEvalN(N, *lcoords, FEIElementGeometryWrapper(this));
            double dV = this->computeEdgeVolumeAround(gp, iEdge);

            if ( boundaryLoad->giveFormulationType() == BoundaryLoad :: BL_EntityFormulation ) {                // Edge load in xi-eta system
                boundaryLoad->computeValueAt(loadValue, tStep, *lcoords, VM_Total);
            } else {  // Edge load in x-y system
                FloatArray gcoords;
                this->interpolation_lin.edgeLocal2global(gcoords, iEdge, *lcoords, FEIElementGeometryWrapper(this));
                boundaryLoad->computeValueAt(loadValue, tStep, gcoords, VM_Total);
            }

            reducedAnswer.add(loadValue.at(1) * dV, N);
        }

        this->interpolation_lin.computeLocalEdgeMapping(mask, iEdge);
        answer.assemble(reducedAnswer, mask);
    }
}

void Tet21Stokes :: computeBoundaryLoadVector(FloatArray &answer, BoundaryLoad *load, int iSurf, CharType type, ValueModeType mode, TimeStep *tStep)
{
    if ( type != ExternalForcesVector ) {
        answer.clear();
        return;
    }

    answer.resize(34);
    answer.zero();

    if ( load->giveType() == TransmissionBC ) { // Neumann boundary conditions (traction)

        int numberOfSurfaceIPs = ( int ) ceil( ( load->giveApproxOrder() + 1. ) / 2. ) * 2; ///@todo Check this.

        GaussIntegrationRule iRule(1, this, 1, 1);
        FloatArray N, t, f(18);

        f.zero();
        iRule.SetUpPointsOnTriangle(numberOfSurfaceIPs, _Unknown);

        for ( GaussPoint *gp: iRule ) {
            const FloatArray &lcoords = gp->giveNaturalCoordinates();

            this->interpolation_quad.surfaceEvalN( N, iSurf, lcoords, FEIElementGeometryWrapper(this) );
            double dA = gp->giveWeight() * this->interpolation_quad.surfaceGiveTransformationJacobian( iSurf, lcoords, FEIElementGeometryWrapper(this) );

            if ( load->giveFormulationType() == Load :: FT_Entity ) { // load in xi-eta system
                load->computeValueAt(t, tStep, lcoords, VM_Total);
            } else { // Edge load in x-y system
                FloatArray gcoords;
                this->interpolation_quad.surfaceLocal2global( gcoords, iSurf, lcoords, FEIElementGeometryWrapper(this) );
                load->computeValueAt(t, tStep, gcoords, VM_Total);
            }

            // Reshape the vector
            for ( int j = 0; j < N.giveSize(); j++ ) {
                f(3 * j + 0) += N(j) * t(0) * dA;
                f(3 * j + 1) += N(j) * t(1) * dA;
                f(3 * j + 2) += N(j) * t(2) * dA;
            }
        }

        answer.assemble(f, this->surf_ordering [ iSurf - 1 ]);
    } else {
        OOFEM_ERROR("Strange boundary condition type");
    }
}