C++ EngngModel::giveDomain方法代码示例

void
StaticFracture :: filterSensitivities(ObjectiveFunction *objFunc)
{
    // Loop through all active elements and compute distance to neighbours

    EngngModel *sp = this->giveSlaveProblem(1);
    Domain *d = sp->giveDomain(1);   
    int numEl = d->giveNumberOfElements(); 
    double radius = 5.0;

    FloatArray oldSensitivities = objFunc->sensitivityList;

    // navie implementation
    FloatArray center0, center, lCoords;
    for (int i = 1; i <= numEl; i++) {
        Element *el = d->giveElement(i);
        lCoords.resize( el->giveSpatialDimension() );
        lCoords.zero(); // center of element
        
        el->computeGlobalCoordinates(center0, lCoords);
        double distSum = 0.0, help = 0.0;
        for (int j = 1; j <= numEl; j++) {
            d->giveElement(j)->computeGlobalCoordinates(center, lCoords);
            double dist = sqrt( center0.distance_square(center) );
            if ( dist <= radius ) {
                distSum += dist;
                help += oldSensitivities.at(j) * dist;
            }
            
        }

        objFunc->sensitivityList.at(i) = help / distSum; // new filtered sensitivity
    }
}

void
SolutionbasedShapeFunction :: giveValueAtPoint(FloatArray &answer, const FloatArray &coords, IntArray &dofIDs, EngngModel &myEngngModel)
{
    answer.resize( dofIDs.giveSize() );

    FloatArray closest, lcoords, values;

    Element *elementAtCoords = myEngngModel.giveDomain(1)->giveSpatialLocalizer()->giveElementClosestToPoint(lcoords, closest, coords, 1);
    if ( elementAtCoords == NULL ) {
        OOFEM_WARNING("Cannot find element closest to point");
        coords.pY();
    }

    EIPrimaryUnknownMapperInterface *em = dynamic_cast< EIPrimaryUnknownMapperInterface * >( elementAtCoords->giveInterface(EIPrimaryUnknownMapperInterfaceType) );

    IntArray eldofids;

    em->EIPrimaryUnknownMI_givePrimaryUnknownVectorDofID(eldofids);
    em->EIPrimaryUnknownMI_computePrimaryUnknownVectorAtLocal(VM_Total, thisTimestep, lcoords, values);

    for ( int i = 1; i <= dofIDs.giveSize(); i++ ) {
        for ( int j = 1; j <= eldofids.giveSize(); j++ ) {
            if ( dofIDs.at(i) == eldofids.at(j) ) {
                answer.at(i) = values.at(j);
                break;
            }
        }
    }
}

void
SolutionbasedShapeFunction :: giveValueAtPoint(FloatArray &answer, const FloatArray &coords, IntArray &dofIDs, EngngModel &myEngngModel)
{
    answer.resize( dofIDs.giveSize() );

    FloatArray closest, lcoords, values;

    Element *elementAtCoords = myEngngModel.giveDomain(1)->giveSpatialLocalizer()->giveElementClosestToPoint(lcoords, closest, coords, 1);
    if ( elementAtCoords == NULL ) {
        OOFEM_WARNING("Cannot find element closest to point");
        coords.pY();
        return;
    }

    IntArray eldofids;

    elementAtCoords->giveElementDofIDMask(eldofids);
    elementAtCoords->computeField(VM_Total, thisTimestep, lcoords, values);

    for ( int i = 1; i <= dofIDs.giveSize(); i++ ) {
        for ( int j = 1; j <= eldofids.giveSize(); j++ ) {
            if ( dofIDs.at(i) == eldofids.at(j) ) {
                answer.at(i) = values.at(j);
                break;
            }
        }
    }
}

void 
StaticFracture :: optimalityCriteria(ObjectiveFunction *objFunc)
{
    FloatArray &x = objFunc->designVarList; 
    FloatArray &dc = objFunc->sensitivityList;
    ///x.printYourself();

    EngngModel *sp = this->giveSlaveProblem(1);
    Domain *d = sp->giveDomain(1);   
    int numEl = d->giveNumberOfElements(); 
    
   
    // Solver paramters
    double l1 = 0; 
    double l2 = 100000; 
    double move = 0.2;
    double averageDensity = 0.0;

    FloatArray xOld = x;
    while (l2-l1 > 1.0e-4) {
        
        double lmid = 0.5*(l2+l1);
        double totalVolume = 0.0;
        double help = 0.0;
        for (int i = 1; i <= numEl; i++) { // should only be el that contribute
            Element *el = d->giveElement(i);
            // xnew = max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid)))));
            double xe = max(0.001, max(xOld.at(i)-move, min(1., min(xOld.at(i) + move, xOld.at(i) * sqrt(-dc.at(i)/lmid) ) )));
            //double xe = max(0.001, max(xOld.at(i)-move, min(1., min(xOld.at(i) + move, xOld.at(i) * sqrt(-objFunc->sensitivityList.at(i)/lmid) ) )));
            //double temp1 = x.at(i) + move;
        //double temp2 = x.at(i) * sqrt(-dc.at(i)/lmid);
        //double temp3 = min(temp1, temp2);
        //temp2 = min(1.0, temp3);
        //temp1 = x.at(i) - move;
        //temp3 = max(temp1, temp2);
        //temp1 = max(0.001,temp3);
        ////designVarList.at(i) = max(0.001, max(x-move, min(1., min(x+move,x.*sqrt(-dCostFunction.at(i)/lmid)))));
            
            x.at(i) = xe;
            //objFunc->designVarList.at(i) = xe;
            double dV = el->computeVolumeAreaOrLength();
            totalVolume += dV;
            help += xe * dV; 
        }
        averageDensity = help / totalVolume; 

        if ( averageDensity - objFunc->constraint > 0 ) {
        //if ( objFunc->designVarList.sum() - objFunc->constraint*20*20 > 0 ) {
            l1 = lmid;
        } else {
            l2 = lmid;
        }
    }

    printf("\n average density %e \n", averageDensity );
}

int
StaticFracture :: instanciateYourself(DataReader *dr, InputRecord *ir, const char *dataOutputFileName, const char *desc)
{

    this->Instanciate_init(dataOutputFileName, this->ndomains);

    fprintf( outputStream, "%s", PRG_HEADER);
    this->startTime = time(NULL);
    fprintf( outputStream, "\nStarting analysis on: %s\n", ctime(& this->startTime) );
    fprintf( outputStream, "%s\n", desc);


    // instanciate receiver
    this->initializeFrom(ir);
    this->instanciateDefaultMetaStep(ir); // after numSteps has been set

    //exportModuleManager->initializeFrom(ir);
    //initModuleManager->initializeFrom(ir);

    int result;
    ir->finish();
    // instanciate slave problems
    result &= this->instanciateSlaveProblems();



    // should loop over slave problems and obj. fncs.
    EngngModel *sp = emodelList->at(1);
    Domain *d = sp->giveDomain(1);   
    int numMat = d->giveNumberOfMaterialModels();


    // initialize Objective function
    MinCompliance *objFunc = dynamic_cast< MinCompliance * >( this->objFuncList[0] ); 

    // Assumed that there is only one design parameter per element
    objFunc->designVarList.resize(numMat);
    objFunc->sensitivityList.resize(numMat);
    for (int i = 1; i <= numMat; i++) {
        objFunc->designVarList.at(i) = objFunc->volFrac;
        objFunc->sensitivityList.at(i) = 0.0;
    }


    return result;
}

void
StaticFracture :: optimize(TimeStep *tStep)
{
    // Main optimization loop

    MinCompliance *objFunc = dynamic_cast< MinCompliance * >( this->objFuncList[0] ); 

    for ( int subProb = 1; subProb <= this->giveNumberOfSlaveProblems(); subProb++ ) {
        
        EngngModel *sp = this->giveSlaveProblem(subProb);
        Domain *d = sp->giveDomain(1);   

        double cost = 0.0; // should lie in obj fnc
        double dce = 0.0;
        for (int i = 1; i <= d->giveNumberOfElements(); i++) {  
            Element *el = d->giveElement(i);
            cost += objFunc->evaluateYourself(el, dce, sp->giveCurrentStep() ); // add cost for each element
        }


        // Filter sensitivities
        this->filterSensitivities(objFunc);

        // Update design variables based on some method. For now use the 'standard' optimality criteria
        this->optimalityCriteria(objFunc);
     

        // Update material parameters
        for (int i = 1; i <= d->giveNumberOfMaterialModels(); i++) {
            DynamicInputRecord ir;
            Material *mat = d->giveMaterial(i);
            mat->giveInputRecord(ir);
            double E0 = 1.0;
            double fac = pow( objFunc->designVarList.at(i), objFunc->penalty);
            ir.setField( E0 * fac, _IFT_IsotropicLinearElasticMaterial_e);
            mat->initializeFrom(&ir);
        }

        printf("\n costfunction %e & sum sensitivity %e & sum x %e \n", cost, objFunc->sensitivityList.sum(), 
            objFunc->designVarList.sum() );

    }
}

void
SolutionbasedShapeFunction :: initializeSurfaceData(modeStruct *mode)
{
    EngngModel *m = mode->myEngngModel;
    double TOL2 = 1e-5;
    IntArray pNodes, mNodes, zNodes;

    Set *mySet = this->domain->giveSet( this->giveSetNumber() );
    IntArray BoundaryList = mySet->giveBoundaryList();

    // First add all nodes to pNodes or nNodes respectively depending on coordinate and normal.
    for ( int i = 0; i < BoundaryList.giveSize() / 2; i++ ) {
        int ElementID = BoundaryList(2 * i);
        int Boundary = BoundaryList(2 * i + 1);

        Element *e = m->giveDomain(1)->giveElement(ElementID);
        FEInterpolation *geoInterpolation = e->giveInterpolation();

        // Check all sides of element
        IntArray bnodes;

#define usePoints 1
#if usePoints == 1
        // Check if all nodes are on the boundary
        geoInterpolation->boundaryGiveNodes(bnodes, Boundary);
        for ( int k = 1; k <= bnodes.giveSize(); k++ ) {
            DofManager *dman = e->giveDofManager( bnodes.at(k) );
            for ( int l = 1; l <= dman->giveCoordinates()->giveSize(); l++ ) {
                if ( fabs( dman->giveCoordinates()->at(l) - maxCoord.at(l) ) < TOL2 ) {
                    pNodes.insertOnce( dman->giveNumber() );
                }
                if ( fabs( dman->giveCoordinates()->at(l) - minCoord.at(l) ) < TOL2 ) {
                    mNodes.insertOnce( dman->giveNumber() );
                }
            }
        }
#else
        // Check normal
        FloatArray lcoords;
        lcoords.resize(2);
        lcoords.at(1) = 0.33333;
        lcoords.at(2) = 0.33333;

        FloatArray normal;
        geoInterpolation->boundaryEvalNormal( normal, j, lcoords, FEIElementGeometryWrapper(e) );
        geoInterpolation->boundaryGiveNodes(bnodes, j);

        printf( "i=%u\tj=%u\t(%f\t%f\t%f)\n", i, j, normal.at(1), normal.at(2), normal.at(3) );
        for ( int k = 1; k <= normal.giveSize(); k++ ) {
            if ( fabs( ( fabs( normal.at(k) ) - 1 ) ) < 1e-4 ) { // Points in x, y or z direction
                addTo = NULL;
                if ( normal.at(k) > 0.5 ) {
                    addTo = & pNodes;
                }
                if ( normal.at(k) < -0.5 ) {
                    addTo = & mNodes;
                }
                if ( addTo != NULL ) {
                    for ( int l = 1; l <= bnodes.giveSize(); l++ ) {
                        bool isSurface = false;
                        DofManager *dman = e->giveDofManager( bnodes.at(l) );
                        dman->giveCoordinates()->printYourself();
                        for ( int m = 1; m <= dman->giveCoordinates()->giveSize(); m++ ) {
                            if ( ( fabs( dman->giveCoordinates()->at(m) - maxCoord.at(m) ) < TOL2 ) || ( fabs( dman->giveCoordinates()->at(m) - minCoord.at(m) ) < TOL2 ) ) {
                                isSurface = true;
                            }
                        }

                        if ( isSurface ) {
                            addTo->insertOnce( e->giveDofManagerNumber( bnodes.at(l) ) );
                        }
                    }
                }
            }
        }
#endif
    }

#if 0
    printf("p=[");
    for ( int i = 1; i < pNodes.giveSize(); i++ ) {
        printf( "%u, ", pNodes.at(i) );
    }
    printf("];\n");
    printf("m=[");
    for ( int i = 1; i < mNodes.giveSize(); i++ ) {
        printf( "%u, ", mNodes.at(i) );
    }
    printf("];\n");
#endif
    //The intersection of pNodes and mNodes constitutes zNodes
    {
        int i = 1, j = 1;
        while ( i <= pNodes.giveSize() ) {
            j = 1;
            while ( j <= mNodes.giveSize() && ( i <= pNodes.giveSize() ) ) {
                //printf("%u == %u?\n", pNodes.at(i), mNodes.at(j));
                if ( pNodes.at(i) == mNodes.at(j) ) {
                    zNodes.insertOnce( pNodes.at(i) );
                    pNodes.erase(i);
//.........这里部分代码省略.........

void
SolutionbasedShapeFunction :: computeBaseFunctionValueAt(FloatArray &answer, FloatArray &coords, IntArray &dofIDs, EngngModel &myEngngModel)
{
    answer.resize( dofIDs.giveSize() );
    answer.zero();

    if ( useConstantBase ) {
        for ( int i = 1; i <= answer.giveSize(); i++ ) {
            answer.at(i) = 1;
        }
        ;
    } else {
        std :: vector< FloatArray * >checkcoords;
        std :: vector< int >permuteIndex;
        int n = 0;

        if ( !isLoaded ) {
            loadProblem();
        }

        myEngngModel.giveDomain(1)->giveSpatialLocalizer()->init(false);

        if ( this->giveDomain()->giveNumberOfSpatialDimensions() == 2 ) {
            coords.resize(2);
        }

        // Determine if current coordinate is at a max or min point and if so, which type (on a surface, edge or a corner?)
        // n is the number of dimensions at which the point is an extremum, permuteIndex tells in which dimension the coordinate is max/min

        int thisMask = 0;

        for ( int i = 1; i <= coords.giveSize(); i++ ) {
            if ( ( fabs( maxCoord.at(i) - coords.at(i) ) < TOL ) || ( fabs( minCoord.at(i) - coords.at(i) ) < TOL ) ) {
                permuteIndex.push_back(i);
                n++;
                //thisMask = thisMask + pow(2.0, i - 1);   // compiler warning on conversion from double to int
                thisMask = thisMask + ( 0x01 << ( i - 1 ) );
            }
        }
        int _s = 0x01 << n;
        for ( int i = 0; i < _s; i++ ) {
            int mask = i, counter = 1;
            FloatArray *newCoord = new(FloatArray) ( coords.giveSize() );
            * newCoord = coords;

            for ( int j = 1; j <= n; j++ ) {
                double d = 0.0; //TOL;
                if ( ( mask & 1 ) == 0 ) { // Max
                    newCoord->at( permuteIndex.at(counter - 1) ) = minCoord.at( permuteIndex.at(counter - 1) ) + d;
                } else { // Min
                    newCoord->at( permuteIndex.at(counter - 1) ) = maxCoord.at( permuteIndex.at(counter - 1) ) - d;
                }
                counter++;
                mask = mask >> 1;
            }
            checkcoords.push_back(newCoord);
        }

        // The followind define allows for use of weakly periodic bc to be copied. This does not comply with the theory but is used to check the validity of the code.
#define USEWPBC 0

#if USEWPBC == 1
        FloatArray *tempCoord = new FloatArray;
        * tempCoord = coords;
        checkcoords.clear();
        checkcoords.push_back(tempCoord);
#endif
        FloatArray values;
        for ( size_t i = 0; i < checkcoords.size(); i++ ) {
            giveValueAtPoint(values, * checkcoords.at(i), dofIDs, myEngngModel);
            //printf("Values at (%f, %f, %f) are [%f, %f, %f]\n", checkcoords.at(i)->at(1), checkcoords.at(i)->at(2), checkcoords.at(i)->at(3), values.at(1), values.at(2), values.at(3));
#if USEWPBC == 1
            for ( int j = 1; j <= values.giveSize(); j++ ) {
                answer.at(j) = values.at(j);
            }
#else
            for ( int j = 1; j <= values.giveSize(); j++ ) {
                answer.at(j) = answer.at(j) + values.at(j) / ( ( double ) pow(2.0, n) );
            }
#endif
            delete( checkcoords.at(i) );
        }
    }
}

void
SolutionbasedShapeFunction :: loadProblem()
{
    for ( int i = 0; i < this->domain->giveNumberOfSpatialDimensions(); i++ ) {
        OOFEM_LOG_INFO("************************** Instanciating microproblem from file %s for dimension %u\n", filename.c_str(), i);

        // Set up and solve problem
        OOFEMTXTDataReader drMicro( filename.c_str() );
        EngngModel *myEngngModel = InstanciateProblem(& drMicro, _processor, 0, NULL, false);
        drMicro.finish();
        myEngngModel->checkProblemConsistency();
        myEngngModel->initMetaStepAttributes( myEngngModel->giveMetaStep(1) );
        thisTimestep = myEngngModel->giveNextStep();
        myEngngModel->init();
        this->setLoads(myEngngModel, i + 1);

        // Check
        for ( int j = 1; j <= myEngngModel->giveDomain(1)->giveNumberOfElements(); j++ ) {
            Element *e = myEngngModel->giveDomain(1)->giveElement(j);
            FloatArray centerCoord;
            int vlockCount = 0;
            centerCoord.resize(3);
            centerCoord.zero();

            for ( int k = 1; k <= e->giveNumberOfDofManagers(); k++ ) {
                DofManager *dman = e->giveDofManager(k);
                centerCoord.add( * dman->giveCoordinates() );
                for ( Dof *dof: *dman ) {
                    if ( dof->giveBcId() != 0 ) {
                        vlockCount++;
                    }
                }
            }
            if ( vlockCount == 30 ) {
                OOFEM_WARNING("Element over-constrained (%u)! Center coordinate: %f, %f, %f\n", e->giveNumber(), centerCoord.at(1) / 10, centerCoord.at(2) / 10, centerCoord.at(3) / 10);
            }
        }

        myEngngModel->solveYourselfAt(thisTimestep);
        isLoaded = true;

        // Set correct export filename
        std :: string originalFilename;
        originalFilename = myEngngModel->giveOutputBaseFileName();
        if ( i == 0 ) {
            originalFilename = originalFilename + "_X";
        }
        if ( i == 1 ) {
            originalFilename = originalFilename + "_Y";
        }
        if ( i == 2 ) {
            originalFilename = originalFilename + "_Z";
        }
        myEngngModel->letOutputBaseFileNameBe(originalFilename + "_1_Base");
        myEngngModel->doStepOutput(thisTimestep);

        modeStruct *mode = new(modeStruct);
        mode->myEngngModel = myEngngModel;

        // Check elements

        // Set unknowns to the mean value of opposite sides of the domain.
        // Loop thru all nodes and compute phi for all degrees of freedom on the boundary. Save phi in a list for later use.

        initializeSurfaceData(mode);
        // Update with factor
        double am = 1.0, ap = 1.0;
        computeCorrectionFactors(* mode, & dofs, & am, & ap);

        OOFEM_LOG_INFO("Correction factors: am=%f, ap=%f\n", am, ap);

        mode->ap = ap;
        mode->am = am;

        updateModelWithFactors(mode);

        //        myEngngModel->letOutputBaseFileNameBe(originalFilename + "_2_Updated");

        modes.push_back(mode);

        OOFEM_LOG_INFO("************************** Microproblem at %p instanciated \n", myEngngModel);
    }
}

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 );
}