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

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 Tr1Darcy :: computeLoadVector(FloatArray &answer, TimeStep *atTime)
{
    // TODO: Implement support for body forces

    FloatArray vec;

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

    // Compute characteristic vector for Neumann boundary conditions.
    int i, load_number, load_id;
    GeneralBoundaryCondition *load;
    bcGeomType ltype;

    int nLoads = boundaryLoadArray.giveSize() / 2;

    for ( i = 1; i <= nLoads; i++ ) {  // For each Neumann boundary condition ....
        load_number = boundaryLoadArray.at(2 * i - 1);
        load_id = boundaryLoadArray.at(2 * i);
        load = ( GeneralBoundaryCondition * ) domain->giveLoad(load_number);
        ltype = load->giveBCGeoType();

        if ( ltype == EdgeLoadBGT ) {
            this->computeEdgeBCSubVectorAt(vec, ( Load * ) load, load_id, atTime);
        }

        answer.add(vec);
    }

    answer.negated();
}

void PrescribedGradientBCWeak :: computeIntForceGPContrib(FloatArray &oContrib_disp, IntArray &oDisp_loc_array, FloatArray &oContrib_trac, IntArray &oTrac_loc_array,TracSegArray &iEl, GaussPoint &iGP, int iDim, TimeStep *tStep, const FloatArray &iBndCoord, const double &iScaleFac, ValueModeType mode, CharType type, const UnknownNumberingScheme &s)
{

    SpatialLocalizer *localizer = domain->giveSpatialLocalizer();

	FloatMatrix contrib;
	assembleTangentGPContributionNew(contrib, iEl, iGP, iScaleFac, iBndCoord);

    // Compute vector of traction unknowns
    FloatArray tracUnknowns;
    iEl.mFirstNode->giveUnknownVector(tracUnknowns, giveTracDofIDs(), mode, tStep);

    iEl.giveTractionLocationArray(oTrac_loc_array, type, s);

    FloatArray dispElLocCoord, closestPoint;
    Element *dispEl = localizer->giveElementClosestToPoint(dispElLocCoord, closestPoint, iBndCoord );

    // Compute vector of displacement unknowns
    FloatArray dispUnknowns;
    int numDMan = dispEl->giveNumberOfDofManagers();
    for(int i = 1; i <= numDMan; i++) {
    	FloatArray nodeUnknowns;
    	DofManager *dMan = dispEl->giveDofManager(i);

    	IntArray dispIDs = giveRegularDispDofIDs();
        if(domain->hasXfemManager()) {
        	XfemManager *xMan = domain->giveXfemManager();
        	dispIDs.followedBy(xMan->giveEnrichedDofIDs(*dMan));
        }

        dMan->giveUnknownVector(nodeUnknowns, dispIDs,mode, tStep);
        dispUnknowns.append(nodeUnknowns);

    }

    dispEl->giveLocationArray(oDisp_loc_array, s);


    oContrib_disp.beTProductOf(contrib, tracUnknowns);
    oContrib_disp.negated();

    oContrib_trac.beProductOf(contrib, dispUnknowns);
    oContrib_trac.negated();
}

void
CBSElement :: computePrescribedTermsII(FloatArray &answer, ValueModeType mode, TimeStep *tStep)
{
    FloatMatrix lhs;
    FloatArray usp;
    this->computePressureLhs(lhs, tStep);
    this->computeVectorOfPressures(mode, tStep, usp);
    answer.beProductOf(lhs, usp);
    answer.negated();
}

void
CBSElement :: computePrescribedTermsI(FloatArray &answer, TimeStep *tStep)
{
    FloatMatrix mass;
    FloatArray usp;
    this->computeConsistentMassMtrx(mass, tStep);
    this->computeVectorOfVelocities(VM_Incremental, tStep, usp);
    answer.beProductOf(mass, usp);
    answer.negated();
}

void
CBSElement :: computePrescribedTermsI(FloatArray &answer, ValueModeType mode, TimeStep *tStep)
{
    FloatMatrix mass;
    FloatArray usp;
    this->computeConsistentMassMtrx(mass, tStep);
    this->computeVectorOf(EID_MomentumBalance, mode, tStep, usp);
    answer.beProductOf(mass, usp);
    answer.negated();
}

void
AnisotropicMassTransferMaterial :: giveFluxVector(FloatArray& answer, GaussPoint *gp, const FloatArray &grad, const FloatArray &field, TimeStep *tStep)
{
    TransportMaterialStatus *ms = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );

    answer.beProductOf(k, grad);
    answer.negated();

    ms->setTempField(field);
    ms->setTempGradient(grad);
    ms->setTempFlux(answer);
}

void
StokesFlowVelocityHomogenization :: applyPressureGradient(const FloatArray &grad)
{
    FloatArray components = grad;
    components.negated();

    ///@todo This should either be a specialized boundary condition so that this routine can be replaced by something else
    for ( auto &bc : this->giveDomain(1)->giveBcs() ) {
        DeadWeight *load = dynamic_cast< DeadWeight* >( bc.get() );
        if ( load ) {
            load->setDeadWeighComponents(components);
            break;
        }
    }
}

void
NonlinearMassTransferMaterial :: giveFluxVector(FloatArray &answer, GaussPoint *gp, const FloatArray &eps, TimeStep *tStep)
{
    AnisotropicMassTransferMaterialStatus *thisMaterialStatus;
    thisMaterialStatus = ( ( AnisotropicMassTransferMaterialStatus * ) this->giveStatus(gp) );

    thisMaterialStatus->setPressureGradient(eps);

    double gradPNorm = eps.computeNorm();

    answer = eps;
    answer.times( 1 + C * pow(gradPNorm, alpha) );
    answer.negated();

    thisMaterialStatus->setSeepageValocity(answer);
}

FloatArray *
J2plasticMaterial :: ComputeStressSpaceHardeningVars(GaussPoint *gp,
                                                     FloatArray *strainSpaceHardeningVariables)
{
    // in full stress strain space

    int i;
    int count = 0, size = this->giveSizeOfFullHardeningVarsVector(), isize, rSize;
    IntArray mask;

    if ( !hasHardening() ) {
        return NULL;
    }

    FloatArray *answer = new FloatArray(size);
    this->giveStressStrainMask( mask, ReducedForm, gp->giveMaterialMode() );
    isize = mask.giveSize();
    rSize = this->giveSizeOfReducedHardeningVarsVector(gp);

    /* kinematic hardening variables are first */
    if ( this->kinematicHardeningFlag ) {
        for ( i = 1; i <= isize; i++ ) {
            // to be consistent with equivalent plastic strain formulation
            // we multiply by (sqrt(2.)*2./3.)
            answer->at( mask.at(i) ) = ( sqrt(2.) * 2. / 3. ) * this->kinematicModuli * strainSpaceHardeningVariables->at(i);
        }

        count = 6;
    }

    if ( this->isotropicHardeningFlag ) {
        answer->at(count + 1) = this->isotropicModuli *
                                strainSpaceHardeningVariables->at(rSize);
    }

    answer->negated();
    return answer;
}

void
J2MPlasticMaterial :: computeStressSpaceHardeningVars(FloatArray &answer, GaussPoint *gp,
                                                      const FloatArray &strainSpaceHardeningVariables)
{
    // in full stress strain space
    int count = 0, size = this->giveSizeOfFullHardeningVarsVector(), isize, rSize;
    IntArray mask;

    if ( !hasHardening() ) {
        answer.clear();
        return;
    }

    answer.resize(size);
    StructuralMaterial :: giveVoigtSymVectorMask( mask, gp->giveMaterialMode() );
    isize = mask.giveSize();
    rSize = this->giveSizeOfReducedHardeningVarsVector(gp);

    /* kinematic hardening variables are first */
    if ( this->kinematicHardeningFlag ) {
        for ( int i = 1; i <= isize; i++ ) {
            // to be consistent with equivalent plastic strain formulation
            // we multiply by (sqrt(2.)*2./3.)
            answer.at( mask.at(i) ) = ( sqrt(2.) * 2. / 3. ) * this->kinematicModuli * strainSpaceHardeningVariables.at(i);
        }

        count = 6;
    }

    if ( this->isotropicHardeningFlag ) {
        answer.at(count + 1) = this->isotropicModuli *
                               strainSpaceHardeningVariables.at(rSize);
    }

    answer.negated();
}

void
LEPlic :: doCellDLS(FloatArray &fvgrad, int ie, bool coord_upd, bool vof_temp_flag)
{
    int i, ineighbr, nneighbr;
    double fvi, fvk, wk, dx, dy;
    bool isBoundaryCell = false;
    LEPlicElementInterface *interface, *ineghbrInterface;
    FloatMatrix lhs(2, 2);
    FloatArray rhs(2), xi(2), xk(2);
    IntArray currCell(1), neighborList;
    ConnectivityTable *contable = domain->giveConnectivityTable();

    if ( ( interface = ( LEPlicElementInterface * ) ( domain->giveElement(ie)->giveInterface(LEPlicElementInterfaceType) ) ) ) {
        if ( vof_temp_flag ) {
            fvi = interface->giveTempVolumeFraction();
        } else {
            fvi = interface->giveVolumeFraction();
        }

        if ( ( fvi > 0. ) && ( fvi <= 1.0 ) ) {
            // potentially boundary cell

            if ( ( fvi > 0. ) && ( fvi < 1.0 ) ) {
                isBoundaryCell = true;
            }

            /* DLS (Differential least square reconstruction)
             *
             * In the DLS method, volume fraction Taylor series expansion of vf (volume fraction)
             * is formed from each reference cell volume fraction vf at element center x(i) to each
             * cell neighbor at point x(k). The sum (vf(i)-vf(k))^2 over all immediate neighbors
             * is then minimized inthe least square sense.
             */
            // get list of neighbours to current cell including current cell
            currCell.at(1) = ie;
            contable->giveElementNeighbourList(neighborList, currCell);
            // loop over neighbors to assemble normal equations
            nneighbr = neighborList.giveSize();
            interface->giveElementCenter(this, xi, coord_upd);
            lhs.zero();
            rhs.zero();
            for ( i = 1; i <= nneighbr; i++ ) {
                ineighbr = neighborList.at(i);
                if ( ineighbr == ie ) {
                    continue;         // skip itself
                }

                if ( ( ineghbrInterface =
                          ( LEPlicElementInterface * ) ( domain->giveElement(ineighbr)->giveInterface(LEPlicElementInterfaceType) ) ) ) {
                    if ( vof_temp_flag ) {
                        fvk = ineghbrInterface->giveTempVolumeFraction();
                    } else {
                        fvk = ineghbrInterface->giveVolumeFraction();
                    }

                    if ( fvk < 1.0 ) {
                        isBoundaryCell = true;
                    }

                    ineghbrInterface->giveElementCenter(this, xk, coord_upd);
                    wk = xk.distance(xi);
                    dx = ( xk.at(1) - xi.at(1) ) / wk;
                    dy = ( xk.at(2) - xi.at(2) ) / wk;
                    lhs.at(1, 1) += dx * dx;
                    lhs.at(1, 2) += dx * dy;
                    lhs.at(2, 2) += dy * dy;

                    rhs.at(1) += ( fvi - fvk ) * dx / wk;
                    rhs.at(2) += ( fvi - fvk ) * dy / wk;
                }
            }

            if ( isBoundaryCell ) {
                // symmetry
                lhs.at(2, 1) = lhs.at(1, 2);

                // solve normal equation for volume fraction gradient
                lhs.solveForRhs(rhs, fvgrad);

                // compute unit normal
                fvgrad.normalize();
                fvgrad.negated();
#ifdef __OOFEG
                /*
                 * EASValsSetLayer(OOFEG_DEBUG_LAYER);
                 * WCRec p[2];
                 * double tx = -fvgrad.at(2), ty=fvgrad.at(1);
                 * p[0].x=xi.at(1)-tx*0.1;
                 * p[0].y=xi.at(2)-ty*0.1;
                 * p[1].x=xi.at(1)+tx*0.1;
                 * p[1].y=xi.at(2)+ty*0.1;
                 * p[0].z = p[1].z = 0.0;
                 * GraphicObj *go = CreateLine3D(p);
                 * EGWithMaskChangeAttributes(LAYER_MASK, go);
                 * EMAddGraphicsToModel(ESIModel(), go);
                 * ESIEventLoop (YES, "Cell DLS finished; Press Ctrl-p to continue");
                 */
#endif
            } else {
                fvgrad.zero();
//.........这里部分代码省略.........

void MixedGradientPressureWeakPeriodic :: assembleVector(FloatArray &answer, TimeStep *tStep,
                                                         CharType type, ValueModeType mode, const UnknownNumberingScheme &s, FloatArray *eNorms)
{
    Set *set = this->giveDomain()->giveSet(this->set);
    const IntArray &boundaries = set->giveBoundaryList();

    IntArray v_loc, t_loc, e_loc;  // For the velocities and stress respectively
    IntArray velocityDofIDs, bNodes;
    this->tractionsdman->giveLocationArray(t_id, t_loc, s);
    this->voldman->giveLocationArray(v_id, e_loc, s);

    if ( type == ExternalForcesVector ) {
        // The external forces have two contributions. on the traction and on dvol.
        double rve_size = this->domainSize();

        if ( e_loc.at(1) ) {
            answer.at( e_loc.at(1) ) -= rve_size * pressure; // Note the negative sign (pressure as opposed to mean stress)
            if ( eNorms ) {
                eNorms->at( v_id.at(1) ) += rve_size * pressure * rve_size * pressure;
            }
        }

        // The second contribution is on the momentumbalance equation; int t . [[ d_dev . x ]] dA = int t . [[ d_dev . x ]] dA
        FloatArray fe;
        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);

            this->integrateTractionDev(fe, e, boundary, this->devGradient);
            fe.negated();

            answer.assemble(fe, t_loc);
            if ( eNorms ) {
                eNorms->assembleSquared(fe, velocityDofIDs);
            }
        }
    } else if ( type == InternalForcesVector ) {
        FloatMatrix Ke_v, Ke_e;
        FloatArray fe_v, fe_t, fe_t2, fe_e(1);
        FloatArray t, e, v;

        // Fetch the current values of internal dofs and their master dof ids;
        this->tractionsdman->giveUnknownVector(t, t_id, mode, tStep);
        this->voldman->giveUnknownVector(e, v_id, mode, tStep);

        // Assemble: -int       t . [[ delta_v ]] dA
        //            int delta_t . [[ e.x - v ]] dA
        //            int       t . [[ x ]]       dA delta_e
        for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
            Element *el = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
            int boundary = boundaries.at(pos * 2);

            // Fetch the element information;
            el->giveInterpolation()->boundaryGiveNodes(bNodes, boundary);
            el->giveBoundaryLocationArray(v_loc, bNodes, this->dofs, s, & velocityDofIDs);
            el->computeBoundaryVectorOf(bNodes, this->dofs, mode, tStep, v);

            // Integrate the tangents;
            this->integrateTractionVelocityTangent(Ke_v, el, boundary);
            this->integrateTractionXTangent(Ke_e, el, boundary);

            // We just use the tangent, less duplicated code
            fe_v.beTProductOf(Ke_v, t);
            fe_v.negated();
            fe_t.beProductOf(Ke_v, v);
            fe_t.negated();
            fe_t2.beProductOf(Ke_e, e);
            fe_t.add(fe_t2);
            fe_e.beTProductOf(Ke_e, t);

            answer.assemble(fe_v, v_loc); // Contributions to delta_v equations
            answer.assemble(fe_t, t_loc); // Contribution to delta_t equations
            answer.assemble(fe_e, e_loc); // Contribution to delta_e equations
            if ( eNorms ) {
                eNorms->assembleSquared(fe_v, velocityDofIDs);
                eNorms->assembleSquared(fe_t, t_id);
                eNorms->assembleSquared(fe_e, v_id);
            }
        }
    }
}

void
NonLinearStatic :: proceedStep(int di, TimeStep *tStep)
{
    //
    // creates system of governing eq's and solves them at given time step
    //
    // first assemble problem at current time step
    //

    if ( initFlag ) {
        //
        // first step  create space for stiffness Matrix
        //
        if ( !stiffnessMatrix ) {
            stiffnessMatrix.reset( classFactory.createSparseMtrx(sparseMtrxType) );
        }

        if ( !stiffnessMatrix ) {
            OOFEM_ERROR("sparse matrix creation failed");
        }

        if ( nonlocalStiffnessFlag ) {
            if ( !stiffnessMatrix->isAsymmetric() ) {
                OOFEM_ERROR("stiffnessMatrix does not support asymmetric storage");
            }
        }

        stiffnessMatrix->buildInternalStructure( this, di, EModelDefaultEquationNumbering() );
    }

#if 0
    if ( ( mstep->giveFirstStepNumber() == tStep->giveNumber() ) ) {
 #ifdef VERBOSE
        OOFEM_LOG_INFO("Resetting load level\n");
 #endif
        if ( mstepCumulateLoadLevelFlag ) {
            cumulatedLoadLevel += loadLevel;
        } else {
            cumulatedLoadLevel = 0.0;
        }
        this->loadLevel = 0.0;
    }
#endif

    if ( loadInitFlag || controlMode == nls_directControl ) {
#ifdef VERBOSE
        OOFEM_LOG_DEBUG("Assembling reference load\n");
#endif
        //
        // assemble the incremental reference load vector
        //
        this->assembleIncrementalReferenceLoadVectors(incrementalLoadVector, incrementalLoadVectorOfPrescribed,
                                                      refLoadInputMode, this->giveDomain(di), tStep);

        loadInitFlag = 0;
    }

    if ( tStep->giveNumber() == 1 ) {
        int neq = this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() );
        totalDisplacement.resize(neq);
        totalDisplacement.zero();
        incrementOfDisplacement.resize(neq);
        incrementOfDisplacement.zero();
    }

    //
    //    ->   BEGINNING OF LOAD (OR DISPLACEMENT) STEP  <-
    //

    //
    // set-up numerical model
    //
    this->giveNumericalMethod( this->giveMetaStep( tStep->giveMetaStepNumber() ) );
    //
    // call numerical model to solve arise problem
    //
#ifdef VERBOSE
    OOFEM_LOG_RELEVANT( "\n\nSolving       [step number %5d.%d, time = %e]\n\n", tStep->giveNumber(), tStep->giveVersion(), tStep->giveIntrinsicTime() );
#endif

    if ( this->initialGuessType == IG_Tangent ) {
#ifdef VERBOSE
        OOFEM_LOG_RELEVANT("Computing initial guess\n");
#endif
        FloatArray extrapolatedForces;
        this->assembleExtrapolatedForces( extrapolatedForces, tStep, TangentStiffnessMatrix, this->giveDomain(di) );
        extrapolatedForces.negated();

        this->updateComponent( tStep, NonLinearLhs, this->giveDomain(di) );
        SparseLinearSystemNM *linSolver = nMethod->giveLinearSolver();
        OOFEM_LOG_RELEVANT("solving for increment\n");
        linSolver->solve(*stiffnessMatrix, extrapolatedForces, incrementOfDisplacement);
        OOFEM_LOG_RELEVANT("initial guess found\n");
        totalDisplacement.add(incrementOfDisplacement);
    } else if ( this->initialGuessType != IG_None ) {
        OOFEM_ERROR("Initial guess type: %d not supported", initialGuessType);
    } else {
        incrementOfDisplacement.zero();
    }

//.........这里部分代码省略.........

void StaticStructural :: solveYourselfAt(TimeStep *tStep)
{
    int neq;
    int di = 1;

    this->field->advanceSolution(tStep);
    this->field->applyBoundaryCondition(tStep); ///@todo Temporary hack, advanceSolution should apply the boundary conditions directly.

    neq = this->giveNumberOfDomainEquations( di, EModelDefaultEquationNumbering() );
    this->field->initialize(VM_Total, tStep, this->solution, EModelDefaultEquationNumbering() ); 

    FloatArray incrementOfSolution(neq), externalForces(neq);

    // Create "stiffness matrix"
    if ( !this->stiffnessMatrix ) {
        this->stiffnessMatrix.reset( classFactory.createSparseMtrx(sparseMtrxType) );
        if ( !this->stiffnessMatrix ) {
            OOFEM_ERROR("Couldn't create requested sparse matrix of type %d", sparseMtrxType);
        }

        this->stiffnessMatrix->buildInternalStructure( this, di, EModelDefaultEquationNumbering() );
    }
    this->internalForces.resize(neq);

    this->giveNumericalMethod( this->giveCurrentMetaStep() );
    this->initMetaStepAttributes( this->giveCurrentMetaStep() );

    if ( this->initialGuessType == IG_Tangent ) {
        OOFEM_LOG_RELEVANT("Computing initial guess\n");
        FloatArray extrapolatedForces;
        this->assembleExtrapolatedForces( extrapolatedForces, tStep, TangentStiffnessMatrix, this->giveDomain(di) );
        extrapolatedForces.negated();

        OOFEM_LOG_RELEVANT("Computing old tangent\n");
        this->updateComponent( tStep, NonLinearLhs, this->giveDomain(di) );
        SparseLinearSystemNM *linSolver = nMethod->giveLinearSolver();
        OOFEM_LOG_RELEVANT("Solving for increment\n");
        linSolver->solve(*stiffnessMatrix, extrapolatedForces, incrementOfSolution);
        OOFEM_LOG_RELEVANT("Initial guess found\n");
        this->solution.add(incrementOfSolution);
    } else if ( this->initialGuessType != IG_None ) {
        OOFEM_ERROR("Initial guess type: %d not supported", initialGuessType);
    } else {
        incrementOfSolution.zero();
    }

    // Build initial/external load
    externalForces.zero();
    this->assembleVector( externalForces, tStep, ExternalForcesVector, VM_Total,
                         EModelDefaultEquationNumbering(), this->giveDomain(1) );
    this->updateSharedDofManagers(externalForces, EModelDefaultEquationNumbering(), LoadExchangeTag);

    if ( this->giveProblemScale() == macroScale ) {
        OOFEM_LOG_INFO("\nStaticStructural :: solveYourselfAt - Solving step %d, metastep %d, (neq = %d)\n", tStep->giveNumber(), tStep->giveMetaStepNumber(), neq);
    }

    double loadLevel;
    int currentIterations;
    NM_Status status = this->nMethod->solve(*this->stiffnessMatrix,
                                            externalForces,
                                            NULL,
                                            this->solution,
                                            incrementOfSolution,
                                            this->internalForces,
                                            this->eNorm,
                                            loadLevel, // Only relevant for incrementalBCLoadVector?
                                            SparseNonLinearSystemNM :: rlm_total,
                                            currentIterations,
                                            tStep);
    if ( !( status & NM_Success ) ) {
        OOFEM_ERROR("No success in solving problem");
    }
}