C++ FEInterpolation::boundaryGiveTransformationJacobian方法代码示例

void
PrescribedMean :: computeDomainSize()
{
    if (domainSize > 0.0) return;


    if (elementEdges) {
        IntArray setList = ((GeneralBoundaryCondition *)this)->giveDomain()->giveSet(set)->giveBoundaryList();

        elements.resize(setList.giveSize() / 2);
        sides.resize(setList.giveSize() / 2);

        for (int i=1; i<=setList.giveSize(); i=i+2) {
            elements.at(i/2+1) = setList.at(i);
            sides.at(i/2+1) = setList.at(i+1);
        }
    } else {
        IntArray setList = ((GeneralBoundaryCondition *)this)->giveDomain()->giveSet(set)->giveElementList();
        elements = setList;
    }

    domainSize = 0.0;

    for ( int i=1; i<=elements.giveSize(); i++ ) {
        int elementID = elements.at(i);
        Element *thisElement = this->giveDomain()->giveElement(elementID);
        FEInterpolation *interpolator = thisElement->giveInterpolation(DofIDItem(dofid));

        IntegrationRule *iRule;

        if (elementEdges) {
            iRule = interpolator->giveBoundaryIntegrationRule(3, sides.at(i));
        } else {
            iRule = interpolator->giveIntegrationRule(3);
        }

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

            double detJ;
            if (elementEdges) {
                detJ = fabs ( interpolator->boundaryGiveTransformationJacobian(sides.at(i), lcoords, FEIElementGeometryWrapper(thisElement)) );
            } else {
                detJ = fabs ( interpolator->giveTransformationJacobian(lcoords, FEIElementGeometryWrapper(thisElement)) );
            }
            domainSize = domainSize + detJ*gp->giveWeight();
        }

        delete iRule;
    }

    printf("%f\n", domainSize);

}

void SurfaceTensionBoundaryCondition :: computeLoadVectorFromElement(FloatArray &answer, Element *e, int side, TimeStep *tStep)
{
    FEInterpolation *fei = e->giveInterpolation();
    if ( !fei ) {
        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);
            }

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

void
SolutionbasedShapeFunction :: computeCorrectionFactors(modeStruct &myMode, IntArray *Dofs, double *am, double *ap)
{
    /*
     * *Compute c0, cp, cm, Bp, Bm, Ap and Am
     */

    double A0p = 0.0, App = 0.0, A0m = 0.0, Amm = 0.0, Bp = 0.0, Bm = 0.0, c0 = 0.0, cp = 0.0, cm = 0.0;

    EngngModel *m = myMode.myEngngModel;
    Set *mySet = m->giveDomain(1)->giveSet(externalSet);

    IntArray BoundaryList = mySet->giveBoundaryList();

    for ( int i = 0; i < BoundaryList.giveSize() / 2; i++ ) {
        int ElementID = BoundaryList(2 * i);
        int Boundary = BoundaryList(2 * i + 1);

        Element *thisElement = m->giveDomain(1)->giveElement(ElementID);
        FEInterpolation *geoInterpolation = thisElement->giveInterpolation();
        IntArray bnodes, zNodes, pNodes, mNodes;
        FloatMatrix nodeValues;

        geoInterpolation->boundaryGiveNodes(bnodes, Boundary);

        nodeValues.resize( this->dofs.giveSize(), bnodes.giveSize() );
        nodeValues.zero();

        // Change to global ID for bnodes and identify the intersection of bnodes and the zero boundary
        splitBoundaryNodeIDs(myMode, * thisElement, bnodes, pNodes, mNodes, zNodes, nodeValues);

        std :: unique_ptr< IntegrationRule >iRule(geoInterpolation->giveBoundaryIntegrationRule(order, Boundary));

        for ( GaussPoint *gp: *iRule ) {
            FloatArray *lcoords = gp->giveCoordinates();
            FloatArray gcoords, normal, N;
            FloatArray Phi;

            double detJ = fabs( geoInterpolation->boundaryGiveTransformationJacobian( Boundary, * lcoords, FEIElementGeometryWrapper(thisElement) ) ) * gp->giveWeight();

            geoInterpolation->boundaryEvalNormal( normal, Boundary, * lcoords, FEIElementGeometryWrapper(thisElement) );
            geoInterpolation->boundaryEvalN( N, Boundary, * lcoords, FEIElementGeometryWrapper(thisElement) );
            geoInterpolation->boundaryLocal2Global( gcoords, Boundary, * lcoords, FEIElementGeometryWrapper(thisElement) );

            FloatArray pPhi, mPhi, zPhi;
            pPhi.resize( Dofs->giveSize() );
            pPhi.zero();
            mPhi.resize( Dofs->giveSize() );
            mPhi.zero();
            zPhi.resize( Dofs->giveSize() );
            zPhi.zero();

            // Build phi (analytical averaging, not projected onto the mesh)
            computeBaseFunctionValueAt(Phi, gcoords, * Dofs, * myMode.myEngngModel);

            // Build zPhi for this DofID
            for ( int l = 1; l <= zNodes.giveSize(); l++ ) {
                int nodeID = zNodes.at(l);
                for ( int m = 1; m <= this->dofs.giveSize(); m++ ) {
                    zPhi.at(m) = zPhi.at(m) + N.at(nodeID) * nodeValues.at(m, nodeID);
                }
            }


            // Build pPhi for this DofID
            for ( int l = 1; l <= pNodes.giveSize(); l++ ) {
                int nodeID = pNodes.at(l);
                for ( int m = 1; m <= this->dofs.giveSize(); m++ ) {
                    pPhi.at(m) = pPhi.at(m) + N.at(nodeID) * nodeValues.at(m, nodeID);
                }
            }

            // Build mPhi for this DofID
            for ( int l = 1; l <= mNodes.giveSize(); l++ ) {
                int nodeID = mNodes.at(l);
                for ( int m = 1; m <= this->dofs.giveSize(); m++ ) {
                    mPhi.at(m) = mPhi.at(m) + N.at(nodeID) * nodeValues.at(m, nodeID);
                }
            }

            c0 = c0 + zPhi.dotProduct(normal, 3) * detJ;
            cp = cp + pPhi.dotProduct(normal, 3) * detJ;
            cm = cm + mPhi.dotProduct(normal, 3) * detJ;

            App = App + pPhi.dotProduct(pPhi, 3) * detJ;
            Amm = Amm + mPhi.dotProduct(mPhi, 3) * detJ;
            A0p = A0p + zPhi.dotProduct(pPhi, 3) * detJ;
            A0m = A0m + zPhi.dotProduct(mPhi, 3) * detJ;

            Bp = Bp + Phi.dotProduct(pPhi, 3) * detJ;
            Bm = Bm + Phi.dotProduct(mPhi, 3) * detJ;
        }
    }

    * am = -( A0m * cp * cp - Bm * cp * cp - A0p * cm * cp + App * c0 * cm + Bp * cm * cp ) / ( App * cm * cm + Amm * cp * cp );
    * ap = -( A0p * cm * cm - Bp * cm * cm - A0m * cm * cp + Amm * c0 * cp + Bm * cm * cp ) / ( App * cm * cm + Amm * cp * cp );
}

void
PrescribedMean :: assemble(SparseMtrx &answer, TimeStep *tStep, CharType type,
                           const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s)
{

    if ( type != TangentStiffnessMatrix && type != StiffnessMatrix ) {
        return;
    }

    computeDomainSize();

    IntArray c_loc, r_loc;
    lambdaDman->giveLocationArray(lambdaIDs, r_loc, r_s);
    lambdaDman->giveLocationArray(lambdaIDs, c_loc, c_s);

    for ( int i=1; i<=elements.giveSize(); i++ ) {
        int elementID = elements.at(i);
        Element *thisElement = this->giveDomain()->giveElement(elementID);
        FEInterpolation *interpolator = thisElement->giveInterpolation(DofIDItem(dofid));

        IntegrationRule *iRule = (elementEdges) ? (interpolator->giveBoundaryIntegrationRule(3, sides.at(i))) :
                                 (interpolator->giveIntegrationRule(3));

        for ( GaussPoint * gp: * iRule ) {
            FloatArray lcoords = gp->giveNaturalCoordinates();
            FloatArray N; //, a;
            FloatMatrix temp, tempT;
            double detJ = 0.0;
            IntArray boundaryNodes, dofids= {(DofIDItem) this->dofid}, r_Sideloc, c_Sideloc;

            if (elementEdges) {
                // Compute boundary integral
                interpolator->boundaryGiveNodes( boundaryNodes, sides.at(i) );
                interpolator->boundaryEvalN(N, sides.at(i), lcoords, FEIElementGeometryWrapper(thisElement));
                detJ = fabs ( interpolator->boundaryGiveTransformationJacobian(sides.at(i), lcoords, FEIElementGeometryWrapper(thisElement)) );
                // Retrieve locations for dofs on boundary
                thisElement->giveBoundaryLocationArray(r_Sideloc, boundaryNodes, dofids, r_s);
                thisElement->giveBoundaryLocationArray(c_Sideloc, boundaryNodes, dofids, c_s);
            } else {
                interpolator->evalN(N, lcoords, FEIElementGeometryWrapper(thisElement));
                detJ = fabs ( interpolator->giveTransformationJacobian(lcoords, FEIElementGeometryWrapper(thisElement) ) );
                IntArray DofIDStemp, rloc, cloc;

                thisElement->giveLocationArray(rloc, r_s, &DofIDStemp);
                thisElement->giveLocationArray(cloc, c_s, &DofIDStemp);

                r_Sideloc.clear();
                c_Sideloc.clear();
                for (int j=1; j<=DofIDStemp.giveSize(); j++) {
                    if (DofIDStemp.at(j)==dofids.at(1)) {
                        r_Sideloc.followedBy({rloc.at(j)});
                        c_Sideloc.followedBy({cloc.at(j)});
                    }
                }
            }

            // delta p part:
            temp = N*detJ*gp->giveWeight()*(1.0/domainSize);
            tempT.beTranspositionOf(temp);

            answer.assemble(r_Sideloc, c_loc, temp);
            answer.assemble(r_loc, c_Sideloc, tempT);
        }

        delete iRule;

    }

}

void
PrescribedMean :: giveInternalForcesVector(FloatArray &answer, TimeStep *tStep,
        CharType type, ValueModeType mode,
        const UnknownNumberingScheme &s, FloatArray *eNorm)
{
    computeDomainSize();

    // Fetch unknowns of this boundary condition
    IntArray lambdaLoc;
    FloatArray lambda;
    lambdaDman->giveUnknownVector(lambda, lambdaIDs, mode, tStep);
    lambdaDman->giveLocationArray(lambdaIDs, lambdaLoc, s);

    for ( int i=1; i<=elements.giveSize(); i++ ) {
        int elementID = elements.at(i);
        Element *thisElement = this->giveDomain()->giveElement(elementID);
        FEInterpolation *interpolator = thisElement->giveInterpolation(DofIDItem(dofid));

        IntegrationRule *iRule = (elementEdges) ? (interpolator->giveBoundaryIntegrationRule(3, sides.at(i))) :
                                 (interpolator->giveIntegrationRule(3));

        for ( GaussPoint * gp: * iRule ) {
            FloatArray lcoords = gp->giveNaturalCoordinates();
            FloatArray a, N, pressureEqns, lambdaEqns;
            IntArray boundaryNodes, dofids= {(DofIDItem) this->dofid}, locationArray;
            double detJ=0.0;

            if (elementEdges) {
                // Compute integral
                interpolator->boundaryGiveNodes( boundaryNodes, sides.at(i) );
                thisElement->computeBoundaryVectorOf(boundaryNodes, dofids, VM_Total, tStep, a);
                interpolator->boundaryEvalN(N, sides.at(i), lcoords, FEIElementGeometryWrapper(thisElement));
                detJ = fabs ( interpolator->boundaryGiveTransformationJacobian(sides.at(i), lcoords, FEIElementGeometryWrapper(thisElement)) );

                // Retrieve locations for dofs with dofids
                thisElement->giveBoundaryLocationArray(locationArray, boundaryNodes, dofids, s);
            } else {
                thisElement->computeVectorOf(dofids, VM_Total, tStep, a);
                interpolator->evalN(N, lcoords, FEIElementGeometryWrapper(thisElement));
                detJ = fabs ( interpolator->giveTransformationJacobian(lcoords, FEIElementGeometryWrapper(thisElement)));

                IntArray DofIDStemp, loc;

                thisElement->giveLocationArray(loc, s, &DofIDStemp);

                locationArray.clear();
                for (int j=1; j<=DofIDStemp.giveSize(); j++) {
                    if (DofIDStemp.at(j)==dofids.at(1)) {
                        locationArray.followedBy({loc.at(j)});
                    }
                }
            }

            // delta p part:
            pressureEqns = N*detJ*gp->giveWeight()*lambda.at(1)*(1.0/domainSize);

            // delta lambda part
            lambdaEqns.resize(1);
            lambdaEqns.at(1) = N.dotProduct(a);
            lambdaEqns.times(detJ*gp->giveWeight()*1.0/domainSize);
            lambdaEqns.at(1) = lambdaEqns.at(1);


            // delta p part
            answer.assemble(pressureEqns, locationArray);

            // delta lambda part
            answer.assemble(lambdaEqns, lambdaLoc);
        }
        delete iRule;
    }

}