• From: Gothenburg, Sweden
  • Registered: 08-09-2009
  • Posts: 1,419

Topic: Homogenization Interface

Me and Carl have concluded roughly what sort of interface we need in order to do homogenization of subscales.
The homogenization does very rarely take part on elements themselves, but typically needs the total stiffness matrix, so one would need a few simple routines on the engineering model, roughly something in the lines of
* Activating homogenization for some certain problem.
* Computing the macro characteristic vector.
* Computing the tangent.

But then I noticed that EngngModel didn't have any giveInterface routine (i guess it doesn't make any sense with any of the current interfaces).

I would add

    /** Returns homogenization interface for given domain */
    virtual HomogenizationInterface *giveHomogenizationInterface() { return NULL; }

to engngm.h

Does that sound acceptable? It didn't make much sense to to ask for a certain domain, so i left that out compared to the other give*** routines.

2 Reply by bp

  • bp
  • Administrator
  • Offline
  • From: Prague, Czech Rep.
  • Registered: 02-03-2007
  • Posts: 1,047

Re: Homogenization Interface

Dear Mikael,

what will the functionality of the interface? To return some characteristic vectors/matrices? In that case I would suggest to integrate these services directly into the class interface, as they will be more likely used by many "interfaces" in the future. On the other hand, if some specific functionality is required (like averaging something, etc) then I would prefer to create specialized derived class for this purpose.  Looking forward to your comments.

Borek

bp's Website

  • From: Gothenburg, Sweden
  • Registered: 08-09-2009
  • Posts: 1,419

Re: Homogenization Interface

I was planning to keep a pure abstract interface.
It is meant to expose a homogenization functionality to material routines.
I was first considering it to be a completely abstract class, leaving it up to engineering models to implement the functionality.

But now I find myself second thinking this design. One could just as well make a different class for each one of these and let engineering models implement whichever they want, since the code for each one (even comparing StressStrainDirichlet with StressStrainNeumann) is very different.
Carl is working on permeability, which has even less in common with what i'm doing when it comes to homogenization.

There is typically no mean value calculation on the elements (but that could also be needed, just very inefficient and complicated in our cases).

As an example, this is what you would do when using dirichlet b.c. on the subscale, in order to calculate the macroscopic stress

void StokesFlow::computeMacroStress(FloatArray &answer, TimeStep *tStep)
{
    // stress = C'.R_c
    int npeq = this->giveNumberOfPrescribedEquations(EID_MomentumBalance);
    FloatArray R_c(npeq), R_ext(npeq);
    R_c.zero(); R_ext.zero();
    this->assembleVectorFromElements( R_c, tStep, EID_MomentumBalance, NodalInternalForcesVector, VM_Total,
                      EModelDefaultPrescribedEquationNumbering(), this->giveDomain(1) );
    // This will also be constant until geometry is updated, I should consider moving this out.
    this->assembleVectorFromElements( R_ext, tStep, EID_MomentumBalance, LoadVector, VM_Total,
                      EModelDefaultPrescribedEquationNumbering(), this->giveDomain(1) );
    R_c.substract(R_ext);

    answer.beTProductOf(*this->C,R_c);
}

and the macroscopic tangent would be some similar matrix manipulation. C is (in 2D) a npeq*3 matrix and contains the coordinates in a strange order (due to Voigt notation). It's quite simple to construct.

The elastic stiffness is simply calculated as

E = C'.(K_cc - K_cf.K_ff^(-1).K_fc).C

where K_xx matrices are from the last iteration of the problem.

So the only thing problem specific is whether EID_MomentumBalance_ConservationEquation or EID_MomentumBalance is used when constructing these matrices.
I could separate the everything non-problem specific and keep these simple matrix operations seperately. The only thing i need from the problem would perhaps be the reaction forces, and all parts of the stiffness matrix.
But some engineering models might be doing some strange things, and is perhaps using some strange numbering schemes.

4 Reply by bp

  • bp
  • Administrator
  • Offline
  • From: Prague, Czech Rep.
  • Registered: 02-03-2007
  • Posts: 1,047

Re: Homogenization Interface

Ok, I see your points. From my point of view, it would be nice if one could implement computeMacroStress method on material level, rather than on problem description level. The reason is to keep the problem representation clean as much as possible. If implemented by some parent class representing material models based on homogenization, then it could be shared by many models. This could be done, provided that the problem can return its characteristic components - like the NodalInternalForcesVector, etc.

It may happen, that some additional functionality will be required by the macro material in addition to methods providing characteristic components. Then the question is whether to use interface concept or make a derived class, that will provide the needed functionality. I would prefer the latter option, as it will allow material to define this specific problem class (ideally as a nested class within itself). The material should be then responsible for maintaining this class. Material also should be responsible for creating this compatible problem representation at micro level.

The interface concept can be used as well. Here I see one main drawback: when several interfaces will come into a play, then there will be a lot of additional services not related to problem solution - this will make the class more difficult not only to follow, but also manage.

What do you think?
Borek

bp's Website

5 Reply by Mikael Öhman (edited by Mikael Öhman 13-06-2010 03:17:35)

  • From: Gothenburg, Sweden
  • Registered: 08-09-2009
  • Posts: 1,419

Re: Homogenization Interface

Edit: I initially replied, but now i have spent some more time thinking this through, and I feel i must change it almost entirely.
If you have read the initial reply already, I apoligize, but I must make changes.

Isolating the homogenization to the material routine, given some interface to obtain characteristic matrices from the problem itself would definitely work for the sort of homogenization i briefly described in my first post. But this is probably the simplest type of homogenization.

But i must make two points to that;
1. Obtaining the matrices is most of the job already. All that is left is to make some simple matrix operations (no more than 10 lines of code, in total for me). Comparing to having to obtain the numerical solver used, load vector, 3 parts of the stiffness matrix (from the last iteration), the coefficient matrix C, i think it would be must simpler to have an interface that returns the homogenized value directly, and a tangent.

2. As you mentioned as well

It may happen, that some additional functionality will be required by the macro material in addition to methods providing characteristic components.

which is what happens if you have anything more complex than my problem.
It typically involves extra equations (mean value of component X should be Y or something similar), which means extra unknowns (typically stresses, traction, or some tensor), and obtaining the tangent might involve having to modify the components in the stiffness matrix.
Some things might also involve integrating over elements.

In short; it can be arbitrary complex.

My old reply, left for historical reasons;

I hadn't thought about isolating it to the material class, but now that you mention it, it does sound like a promising solution.

However, I have some concerns.
Right now the engineering model have very large freedom in how to solve the problem. It could use different equation numbering schemes, it chooses what sort of solver should be used and so on. Just obtaining the stiffness matrices doesn't really give sufficient information for the material routine.

Looking at how I calculate the macro tangent;

E = C'.(K_cc - K_cf.K_ff^(-1).K_fc).C

Here, K_fc.C = neq*3 matrix, and i need to solve for K_ff 3 times. Thus, i also need the nMethod from the problem so solve these three linear problems.

Also, constructing C. Right now, I'm looping over all dofs, looking for which are prescribed, and what dof type is it (D_u/V_u, D_v/V_v, D_w/V_w), and it's prescribed equation number determines the row. It's quite simple, but I need access to the domain, and the nodes to do so.

And my homogenization case is probably the cleanest and simplest type there is.  With different boundary conditions this could become very complex.

Typically when using neumann boundary conditions, this become alot more messy, having to introduce new types of unknowns (ranging from unknown macroscopic traction, or gradients of some macroscale tensor).
For example, some types of boundary condition are applied as
t = g.x
where t = traction, x = coordinate, and g is unknown, to satisfy some extra condition which typically involves from integral over the area.
If I want to use neumann boundary condition for stress/strain, I would get 3 extra unkowns (macroscopic stresses, sigma_11, sigma_22, sigma_12), and some extra equations to satisfy.

It would require so much information from the engineering model, it somewhat loses it's purpose to do it outside.

6 Reply by Mikael Öhman (edited by Mikael Öhman 14-06-2010 22:43:52)

  • From: Gothenburg, Sweden
  • Registered: 08-09-2009
  • Posts: 1,419

Re: Homogenization Interface

Sorry for being a bit talkative, but I've spent some more time trying to design something general, clean and useful, and now I feel a have some more insight. I shall try to give a summarized answer to each part.

bp wrote:

Ok, I see your points. From my point of view, it would be nice if one could implement computeMacroStress method on material level, rather than on problem description level. The reason is to keep the problem representation clean as much as possible. If implemented by some parent class representing material models based on homogenization, then it could be shared by many models. This could be done, provided that the problem can return its characteristic components - like the NodalInternalForcesVector, etc.

Other than constructing the necessary matrices there would only be about 7 short and simple lines of code left for me.
I would have to write more code to return every matrix (and numerical method to solve) then to actually just calculate the stress and tangent directly. Most code lies in constructing the coefficient matrix C, which i could write in a support routine in the interface for the engineering model.

Only this type of homogenization would need these matrices, vectors and objects. A different method would have to do vastly different things. Perhaps iterate over solveYourselfAt(tStep);

It may happen, that some additional functionality will be required by the macro material in addition to methods providing characteristic components.

I'd just like to point out that this is no "what if". This will definitely happen, and not in exotic boundary cases (the 3-4 kinds of permeability homogenization that Carl is working with for example).


Then the question is whether to use interface concept or make a derived class, that will provide the needed functionality. I would prefer the latter option, as it will allow material to define this specific problem class (ideally as a nested class within itself). The material should be then responsible for maintaining this class. Material also should be responsible for creating this compatible problem representation at micro level.

The interface concept can be used as well. Here I see one main drawback: when several interfaces will come into a play, then there will be a lot of additional services not related to problem solution - this will make the class more difficult not only to follow, but also manage.

What do you think?
Borek

In a the derived class, for more complicated homogenization, there will probably be a need to change the way the problem is solved. So, I'm afraid it would lead to having to replace most of the code in the parent class.

The way i see it, you can't easily separate problem solving from homogenization. It's an integral part of solving, as in some cases, it would radically change the problem itself.

Stress/Strain with Dirichlet boundary conditions was probably not a good representative example, although it is probably the most common one.

Yes, it probably make the engineering model more complex for those who do not need homogenization, but isn't the same true for every optional feature?


I would still suggest this interface

/**
 * @file HomogenizationInterface.h
 * @since 2.0
 * @date 2010-06-08
 * @author Mikael Öhman
 */

#ifndef HOMOGENIZATIONINTERFACE_H
#define HOMOGENIZATIONINTERFACE_H

namespace oofem {

enum HomogenizationType {
    HT_None = 0,
    HT_StressStrainDirichlet = 1, /// Normal dirichlet BC for small deformation (structural problems), or velocity (flow problems).
    HT_StressStrainNeumann = 2, /// Normal dirichlet BC for small deformation (structural problems), or velocity (flow problems).
    HT_PermeabilityConstantVelocity = 3, /// For permeability problem, @f$v_i = K_{ij}p_{,j}@f$ using constant velocity on the boundary.
    HT_PermeabilityConstantTraction = 4, /// For permeability problem, @f$v_i = K_{ij}p_{,j}@f$ using constant traction on the boundary.
    HT_PermeabilityLinearTraction = 5 /// For permeability problem, @f$v_i = K_{ij}p_{,j}@f$ using linearly varying traction on the boundary.
};

/**
 * The homogenization interface lets engineering models implement some types of homogenization over the domain.
 * The engineering models can then be used by constitutive drivers to solve micro problems for each integration point.
 */
class HomogenizationInterface : public Interface
{

protected:
    HomogenizationType ht;

public:

    // Perhaps seperate this, and use "solveFor"
    //virtual void solveFor(const FloatArray &input, TimeStep *tStep);

    // Should these also have an input for homogenizationtype? I don't think one would ever need to have more than one at the same time.
    /** Computes the characteristic vector for the active RVE problem.
     * This method should set the boundary condition given the input, solve the micro problem and perform homogenization.
     * @param answer The characteristic vector.
     * @param input Input to set the boundary condition.
     * @param tStep The timestep to solve for.
     * @return True for succes, otherwise false.
     */
    virtual bool computeMacroCharacteristicVector(FloatArray &answer, const FloatArray &input, TimeStep *tStep) = 0;

    /** Computes the tangent for the last calculated vector.
     * computeMacroCharacteristicVector should be called before calling this.
     * @param answer The tangent.
     * @param tStep The timestep to solve for.
     * @return True for success, otherwise false.
     */
    virtual bool computeMacroTangent(FloatMatrix &answer, TimeStep *tStep) = 0;

    /** Prepares to use engineering model as an RVE.
     * Implementation should show error if the type is not supported.
     * @param ht Sets which type of problem/boundary condition should be used for the simulation.
     * @return bool True if type is supported, otherwise false.
     */
    virtual bool activateHomogenizationMode(HomogenizationType ht) = 0;

    virtual const char *giveClassName() const  { return "HomogenizationInterface"; }
};

} // end namespace oofem
#endif // HOMOGENIZATIONINTERFACE_H

There are several ways to handle common operations, but I'm not sure that are so many common things within homogenization.

I could add this add some static support methods so one doesn't need to implement it again.

FloatMatrix* constructCoefficientMatrixDirichlet(Domain *domain, FloatArray &cCoords);
{
    int nNodes = domain->giveNumberOfDofManagers();

    if (domain->giveNumberOfElements() == 0) {
        this->C = new FloatMatrix(0,0);
        return;
    }
    Element *e1 = domain->giveElement(1);
    int dim = e1->giveSpatialDimension();

    int dofType[3] = {0,0,0};
    domainType d = domain->giveDomainType();
    if (d == _2dIncompressibleFlow || d == _3dIncompressibleFlow) {
        dofType[0] = V_u; dofType[1] = V_v; dofType[2] = V_w;
    }
    else if (d == _3dMode || d == _2dPlaneStressMode || d == _PlaneStrainMode || d == _2dTrussMode) {
        dofType[0] = D_u; dofType[1] = D_v; dofType[2] = D_w;
    }
    else
        OOFEM_ERROR("Unsupported domain type to construct C matrix");

    int npeq = this->giveNumberOfPrescribedEquations(EID_MomentumBalance);
    C = new FloatMatrix(npeq, dim*(dim+1)/2);
    C->zero();

    PrescribedTensor *pt = dynamic_cast<PrescribedTensor*>(domain->giveBc(1));
    if (!pt) {
        OOFEM_WARNING("Incorrect boundary condition, should be prescribed tensor");
        return false;
    }

    double xbar = cCoords.at(1), ybar = cCoords.at(2), zbar = cCoords.at(3);
    for (int i = 1; i <= nNodes; i++) {
        Node *n = domain->giveNode(i);
        FloatArray *coords = n->giveCoordinates();
        Dof *d1 = n->giveDofWithID(dofType[0]);
        Dof *d2 = n->giveDofWithID(dofType[1]);
        int k1 = d1->__givePrescribedEquationNumber();
        int k2 = d2->__givePrescribedEquationNumber();
        if (dim == 2) {
            if (k1) {
                C->at(k1,1) = coords->at(1)-xbar;
                C->at(k1,3) = coords->at(2)-ybar;
            }
            if (k2) {
                C->at(k2,2) = coords->at(2)-ybar;
                C->at(k2,3) = coords->at(1)-xbar;
            }
        }
        else { // dim == 3
            OOFEM_ERROR("COMPLETELY UNTESTED!");
            Dof *d3 = n->giveDofWithID(dofType[2]);
            int k3 = d3->__givePrescribedEquationNumber();

            if (k1) {
                C->at(k1,1) = coords->at(1)-xbar;
                C->at(k1,4) = coords->at(2)-ybar;
                C->at(k1,5) = coords->at(3)-zbar;
            }
            if (k2) {
                C->at(k2,2) = coords->at(2)-ybar;
                C->at(k2,4) = coords->at(1)-xbar;
                C->at(k2,6) = coords->at(3)-zbar;
            }
            if (k3) {
                C->at(k3,3) = coords->at(3)-zbar;
                C->at(k3,5) = coords->at(1)-xbar;
                C->at(k3,6) = coords->at(2)-ybar;
            }
        }
    }
    return C;
}

Of course, I have made sure that everything related to homogenization is neatly isolated into functions as much as possible.

I can't see any downsides of doing it this way.

Edit: Had to edit and add text again, the forum tends to cut my replies short.

7 Reply by bp

  • bp
  • Administrator
  • Offline
  • From: Prague, Czech Rep.
  • Registered: 02-03-2007
  • Posts: 1,047

Re: Homogenization Interface

Dear Mikael,

I am sorry for late reply, but I am quite busy at the end of this term.
I still would like to advocate the use of derived classes intended to support homogenization rather than modifying the original plain solver.

As we both agreed, the support for any extension (like homogenization) would require to change the implementation of the corresponding EngngModel. If one thinks in terms of supporting a single extension, than one can perhaps modify the problem implementation. But once you admit, that several extension may exist, then you may quickly end up with very unreadable code as every extension may require to change solution sequence, etc.
From my point of view, the original code, solving the plain problem should be preserved - as it cleanly and efficiently solves the problem.
And I believe, that using derived classes to implement extensions can bring more benefits, even though there will probably be a need to change the way the problem is solved. You have written that it would lead to having to replace most of the code in the parent class. Probably yes, but we can break up the implementation of parent class into smaller methods with the hope of reuse when implementing extension in derived class.
For me it is better to have rather two classes clearly implementing the required stuff, rather than having a single messy one.

Borek

bp's Website

  • From: Gothenburg, Sweden
  • Registered: 08-09-2009
  • Posts: 1,419

Re: Homogenization Interface

Dear Borek,
I was talking to Carl, and he has managed to isolate his homogenization quite well.
I think several derived classes could be the way to go after all.

I'm preparing some patches, but I will leave for a while on thursday, and i don't want to rush you with unpolished patches, so i'll probably won't send you anything until 2 weeks from now.