Interface your own physical solver

This section aims to adress concerns for a developer and does therefore require a deeper understanding of the underlying implementation details.

How to interface your own physical solver to run a physical analysis

The aim of this section is to explain how to run a physical analysis using your own solver. Some aspects related to this question are already partially covered in the Physical Solvers documentation, it is strongly recommended to read it conjointly.

Principles

As of now, there are four interfaces with physical solvers available. All the current interfaces inherit from the same base class and involve the modules OpenPisco.PhysicalSolvers and OpenPisco.ExternalTools. In particular, they provide support for preprocessing and methods to handle the concept of auxiliary fields. In order to set these ideas down, here is an implementation example

from OpenPisco.PhysicalSolvers.SolverBase import SolverBase
import OpenPisco.PhysicalSolvers.FieldsNames as FN

class MyAwesomeSolver(SolverBase):
    def __init__(self):
        super(MyAwesomeSolver,self).__init__()

        #Optional
        self.auxiliaryScalarGeneration[FN.MyScalar] = False
        self.auxiliaryFieldGeneration[FN.MyAuxiliaryField1][FN.Nodes] = False
        self.auxiliaryFieldGeneration[FN.MyAuxiliaryField2][FN.Nodes] = False

    def SolveByLevelSet(self, levelset):
        #Solve physical problem
        return RETURN_SUCCESS

    def GetNodalSolution(self):
        #Retrieve solution at nodes
        raise nodalSolution

    def GetAuxiliaryField(self,name,on):
        #Retrieve auxiliary field
        return auxiliaryFieldOn

    def GetAuxiliaryScalar(self,name):
        #Retrieve auxiliary scalar
        return auxiliaryScalar

In order to interface a new physical solver, a developer should redefine the behaviour of few methods.

Table 5 Methods to be redefined to interface a new physical solver

Method

Notes

OpenPisco.PhysicalSolvers.SolverBase.SolveByLevelSet()

run the analysis

OpenPisco.PhysicalSolvers.SolverBase.GetNodalSolution()

retrieve the PDE solution at the nodes

OpenPisco.PhysicalSolvers.SolverBase.GetAuxiliaryField()

(Optional) Allows a computed auxiliary field to be retrieved

OpenPisco.PhysicalSolvers.SolverBase.GetAuxiliaryScalar()

(Optional) Allows a computed auxiliary scalar to be retrieved

Note

The GetAuxiliary methods implementation depends on your business logic. For instance, for the Code_Aster interface, auxiliary fields are written in the output file produced by Code_Aster. In that case, you have to fulfill accordingly the dictionary structure in the derived class constructor to specify which auxiliary quantities your solver can compute. The auxiliary fields are defined in the OpenPisco.PhysicalSolvers.FieldsNames.

In practice, the solver interface implementation offer therefore some leeway to the user, depending on whether there is a python API:

  • Existing python API: Such is the case, for example, for the solver provided within the Muscat library. This is the most straightforward case to interface as there is no need to build a wrapper nor writing files on disk to use this solver. Not that the last point is not necessarily always true for all solver with a Python API, it turns out the mesh object used through the platform is also a Muscat object.

  • No existing python API: the tool is available by means of an executable. The general procedure in this case is the following:

    1. Writing the input FE model: export the mesh and fields to be used as an input in a suitable format for the solver

    2. Business logic definition: Scripts relying on the solver Domain Specific Langage (DSL) are generated, either from scratch or based on preexisting template available on the platform. If some auxiliary quantities were defined in the problem and need to be computed, it is to be handled.

    3. Run the problem: for instance, building the command line and execute it

    4. Post processing: retrieve the physical problem solution and auxiliary fields, if any

It is recommanded to have a look at the current solver implementations in the dedicated module OpenPisco.PhysicalSolvers:

Compatibility with OpenPisco and OpenPiscoCL applications

The only action you have left is to define how the keywords used within the input file for the applications, relying on OpenPisco DSL, match your own implementation. We shall consider Code_Aster interfacing to illustrate that. Consider the following code from the tutorial Topology optimization of a cantilever beam using unstructured conformal level sets

<PhysicalProblems>
    <GeneralAster type="static_elastic" id="1" >
        <Material  E="210.e9" Nu="0.3" />
        <Dirichlet eTag="eTag1" dofs="0 1 2" value="0.0"/>
        <Dirichlet eTag="eTag3" dofs="1" value="0.0"/>
        <Force nTag="nTag4" value="0. 0. -10000" />
    </GeneralAster>
</PhysicalProblems>

Internally, this code becomes a mere dictionary given to a specific constructor function to initialize properly the physical problem (type of analysis, material properties, boundary conditions…). The function responsible to ensure that for Code_Aster is OpenPisco.PhysicalSolvers.GeneralAsterPhysicalSolver.CreateAsterPhysicalProblemMeca() and return an instance of the suitable physical problem, properly initialized. As this part is completely decoupled from the rest of the input block of code above, with the exception of the id, there is no restriction whatsoever regarding the keywords you use for your own solver.

Now the last step is simply to add to the physical factory the mapping between the solver name, the associated class and associated constructor function. For instance

from OpenPisco.PhysicalSolvers.PhysicalSolverFactory import RegisterClass
RegisterClass("GeneralAster", AsterStatic,CreateAsterPhysicalProblemMeca)

where:

  • “GeneralAster” is the name used in the input block of code above, for the applications

  • AsterStatic is the actual class OpenPisco.PhysicalSolvers.AsterStatic associated to the “static_elastic” analysis

  • CreateAsterPhysicalProblemMeca, mentioned above, is the general constructor-like function for all physical analysis using Code_Aster. Note that its implementation depends on the underlying solver interface implementation.

How to implement your own criterion for topology optimization purposes

We strongly recommanded to read the documentation related to the criteria conjointly.

Principles

All the criteria are available in the module OpenPisco.Optim.Criteria. For practical purposes, there are two main types of criteria considered:

The main difference between them is the endowed physical problem(s) OpenPisco.Optim.Criteria.Criteria.PhysicalCriteriaBase.problem within a physical criterion, which must be run in order to compute the value of the criteria. Typically, its value arise from a postprocessing of the underlying solution to a PDE systems, meaning that such an attribute is an instance of a physical solver supported by the platform. For the sake of completeness, note that a geometrical criterion could also involve the resolution of a PDE problem, as long as it is related to the actual shape of the structure to be optimized.

In order to illustrate that, here is an implementation example for a typical physical criterion

from OpenPisco.Optim.Criteria.CriteriaFactory import RegisterClass as RegisterCriteriaClass

class MyAwesomeCriteria(PhysicalCriteriaBase):
    def __init__(self, other=None):
        super(MyAwesomeCriteria).__init__(other)

    def UpdateValues(self,levelSet):
        #Compute criteria value
        self.f_val=criteriaValue
        #Compute criteria sensitivity
        self.fSensitivity_val=sensitivityValue
        return RETURN_SUCCESS

    def GetNumberOfSolutions(self):
        #Provide numberOfSolution
        return numberOfSolution

    def GetSolution(self,i):
        if i >= self.GetNumberOfSolutions():
            raise(Exception("i out of bounds "))

        #Retrieve solution i
        return ithSolution

    def GetSolutionName(self,i):
        if i >= self.GetNumberOfSolutions():
           raise(Exception("i out of bounds "))

        #Retrieve ith solution name
        return ithSolutionName

Before proceeding, the concept of solutions within a criterion should be explained. A given criterion can involve several fields of interest, including the underlying solutions to the PDE system hidden in the OpenPisco.Optim.Criteria.Criteria.PhysicalCriteriaBase.problem which should the first field of interest to consider in the general case. Throught the optimization process, it can be interesting to export such fields in order to observe how the evolution on the shape impacts them; this is why we provide this possibility for each criteria.

The constraints and responsabilities for the developer are as follows:

  • Redefinition of the behaviour of OpenPisco.Optim.Criteria.Criteria.CriteriaBase.UpdateValues(): it is the main method. This method is the one actually called by the optimization problem (OpenPisco.Optim.Problems.OptimProblemConcrete.UpdateValues()) to update the criteria value/sensitivities. Note that it must return the criterion state at the end of the computation, here RETURN_SUCCESS, if it suceeds.

  • Redefine the behaviour of OpenPisco.Optim.Criteria.Criteria.CriteriaBase.GetNumberOfSolutions(): return the number of field of interest you would like to export

  • Redefine the behaviour of OpenPisco.Optim.Criteria.Criteria.CriteriaBase.GetSolution(): return a field of interest (array) you would like to export for a given index

  • Redefine the behaviour of OpenPisco.Optim.Criteria.Criteria.CriteriaBase.GetSolutionName():return the name of a field of interest (array) you would like to export for a given index

In order to add a new optimization criterion, a developer should redefine the behaviour of few methods.

Table 6 Methods to be redefined to add a new optimization criterion

Method

Notes

OpenPisco.Optim.Criteria.Criteria.CriteriaBase.UpdateValues()

Update the criteria value/sensitivities

OpenPisco.Optim.Criteria.Criteria.CriteriaBase.GetNumberOfSolutions()

return the number of field of interest to export

OpenPisco.Optim.Criteria.Criteria.CriteriaBase.GetSolution()

return a field of interest (array) to export

OpenPisco.Optim.Criteria.Criteria.CriteriaBase.GetSolutionName()

return the name of a field of interest (array) to export

Note

The number of solutions provided should be consistent between these three methods (same number of names, same number of fields and correct number of solution provided).

Compatibility with OpenPisco and OpenPiscoCL applications

Let us consider an example from the tutorial Topology optimization of a cantilever beam using unstructured conformal level sets with the compliance criterion.

<OptimProblems>
    <OptimProblem type="TopoGeneric" id="1"
        OnZone="8" Dirichlet="2" Neumann="5"
        useLevelset="1"
    printMeshQualityInfo = "True"
        output="3"
        outputEveryLevelsetUpdate = "False">
      <Objective  type="Volume"  name="Volume"/>
      <Constraint type="Compliance"  name="Compliance" UpperBound="2.5*4.69434000e-01"   useProblem="1" />
   </OptimProblem>
</OptimProblems>

In order to enable the use of a criterion in the application, in the OpenPisco.Optim.Criteria.PhyCriteria module, we add the mapping between the criterion name and the associated class to the criteria factory.

from OpenPisco.Optim.Criteria.CriteriaFactory import RegisterClass as RegisterCriteriaClass
RegisterCriteriaClass("Compliance",TopoCriteriaCompliance)

where:

  • “Compliance” is the name used in the input block of code above, for the applications

  • TopoCriteriaCompliance is the class OpenPisco.Optim.Criteria.PhyMecaCriteria.TopoCriteriaCompliance associated to the compliance criteria

Note only the simplest configuration is covered above; it is implicitly assumed that the name of the criterion is enough to build the associated criterion instance. Let us consider another criterion in order to see how to handle more complex cases

<OptimProblems>
    <OptimProblem type="TopoGeneric" id="1"
        OnZone="8" Dirichlet="2" Neumann="5"
        useLevelset="1"
    printMeshQualityInfo = "True"
        output="3"
        outputEveryLevelsetUpdate = "False">
      <Objective  type="Volume"  name="Volume"/>
      <Constraint type="NodalTargetDisp"  name="NodalTargetDisp" TargetValue="1.e-3" useProblem="1" dir="0 0 1" u0="0.5" nTag="NT3"/>
   </OptimProblem>
</OptimProblems>

Just like the physical solver, such a block of code becomes internally a mere dictionary. The values within it have to be passed to the criterion instance for its proper initialization. In order to do that, unlike the previous case, we need the associated constructor-like function. For instance

from OpenPisco.Optim.Criteria.CriteriaFactory import RegisterClass as RegisterCriteriaClass
RegisterCriteriaClass("NodalTargetDisp", TopoCriteriaNodalTargetDisp,CreateTopoCriteriaNodalTargetDisp)

where: