ENigMA - Extended Numerical Multiphysics Analysis

Build Status	License

ENigMA - Extended Numerical Multiphysics Analysis

ENigMA is an object-oriented C++ template library which aim is to provide multi-physics simulation in a multi-domain environment. The code implements several numerical methods such as Finite Volume Methods (FVM), Finite Difference Methods (FDM), Finite Element Methods (FEM), Boundary Element Methods (BEM), Smoothed Particle Hydrodynamics (SPH), etc. for numerical approximation of Partial Differential Equations (PDE) in each domain. It also provides classes for robust mesh generation (triangular, block, constrained tetrahedral, etc), intersection and clipping operations and implements R-tree, octree and hashgrid methods for spatial searching. It can be used for three-dimensional flow, thermal and structural analysis.

It was developed to be cross-platform using STL, Eigen (for vectors and matrices) and exprtk (for math expressions). The SWIG tool is used to expose ENigMA's classes and methods to other languages such as Python, C#, etc. It also uses Gtest for unit testing (> 160 unit tests), CMake is used for cross-platform building and Git is used for source code management.

ENigMA implements several operators such as divergence, laplacian, gradient, etc. so, for example, a heat equation such as:

can be represented in ENigMA as:

CPdeEquation<double> aPdeEquation{1.0 / dt * ddt<double>(u) - alpha * laplacian<double>(u) = 0};

The argument passed to the CPdeEquation accepts an object of type CSleSystem which is a linear system of equations. The solution of the PDE is achieved by simply calling solve method:

aPdeEquation.solve(u);

Repository

https://github.com/bjaraujo/ENigMA

Structure

ENigMA has the following namespaces:

analytical: analytical solutions
geometry: geometry, intersection, clipping, spatial search (R-tree, octree, hashgrid), etc.
integration: numerical integration
bem: Boundary Element Method (BEM)
fdm: Finite Difference Method (FDM)
fem: Finite Element Method (FEM)
fvm: Finite Volume Method (FVM)
sph: Smoothed Particle Hydrodynamics (SPH)
pde: Partial Differential Equations
sle: System of Linear Equations
mesh: mesh generation (block, triangular, tetrahedral, hexahedral, etc.)
post: post-processing
stl: STL file processing

Capabilities

Mesh generation (2D/3D),
Structural mechanics (linear elasticity),
Electro-magnetics (Laplace equation),
Fluid flow and heat transfer (Navier-Stokes & Laplace equation),
Finance (Black-Scholes equation).

Dependencies

Install

Configure correctly the environment variables (EIGEN_DIR, EXPRTK_DIR, RTREE_DIR, etc.) and run CMake.

git clone https://github.com/bjaraujo/ENigMA.git
cd ENigMA
cmake-gui&
setup source/build directory
set paths to Exprtk (dependencies folder), RTree (dependencies folder)
configure/generate
make

Usage

Check the tests and examples folders,
Visit the Wiki,
Run demo (releases).

Quick Start

Install python 3 32/64bit (if you don't have already),
Download the latest release,
Install using pip:

pip3 install ENigMApy-VERSION-py3-none-DISTRO.whl

Create file named "Example1.py"
Go to wiki and copy example "Creating a mesh using python"
Run command "python Example1.py"

Examples

These examples show the mesh generation capability of the ENigMA library. The pythonocc library is used to build the CAD geometries and the ENigMA python wrapper to perform surface mesh generation. These examples use Gmsh for post-processing.

To reproduce these examples:

Download miniconda (Python 3.9.12 64bit): https://conda.io/miniconda.html
Install pythonocc:

conda install -c conda-forge pythonocc-core=7.5.1

Install ENigMA 0.4.7:

pip install ENigMApy_0.4.7.0-py3.9-none-win_amd64.whl

Download the examples: ENigMApy_0.4.7.0-py3-occ-examples.zip

Note: If you are using Windows and you get the following error: "ImportError: DLL load failed: The specified module could not be found", you might need to download and install the Microsoft Visual C++ Redistributable for Visual Studio 2015, 2017 and 2019.

Mesh Generation

A Cylinder

Code

from OCC.Core.BRepPrimAPI import BRepPrimAPI_MakeCylinder

import ENigMAocc

cylinder = BRepPrimAPI_MakeCylinder(3.0, 10.0).Shape() 

mesh = ENigMAocc.meshShape(cylinder, 0.5, 1E-3)
ENigMAocc.saveMeshFile(mesh, "occ_01.msh")

A Torus

Code

from OCC.Core.BRepPrimAPI import BRepPrimAPI_MakeTorus

import ENigMAocc
    
sphere = BRepPrimAPI_MakeTorus(4.0, 1.5).Shape() 

mesh = ENigMAocc.meshShape(sphere, 0.5, 1E-3)
ENigMAocc.saveMeshFile(mesh, "occ_02.msh")

A Cut-out Sphere

Code

from OCC.Core.gp import gp_Pnt2d, gp_Pnt, gp_Vec, gp_Trsf
from OCC.Core.BRepPrimAPI import BRepPrimAPI_MakeSphere
from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_Transform
from OCC.Core.BRepAlgoAPI import BRepAlgoAPI_Cut

import ENigMAocc
        
sphere1 = BRepPrimAPI_MakeSphere(3.0).Shape() 
sphere2 = BRepPrimAPI_MakeSphere(6.0).Shape() 

move = gp_Trsf()
move.SetTranslation(gp_Vec(4.0, 0.0, 0.0))
sphere2 = BRepBuilderAPI_Transform(sphere2, move, False).Shape()

shape = BRepAlgoAPI_Cut(sphere1, sphere2).Shape()

mesh = ENigMAocc.meshShape(shape, 0.5, 1E-3)
ENigMAocc.saveMeshFile(mesh, "occ_03.msh")

A Cut-out Shape

Code

from OCC.Core.gp import gp_Vec, gp_Trsf
from OCC.Core.BRepPrimAPI import BRepPrimAPI_MakeBox, BRepPrimAPI_MakeCylinder
from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_Transform
from OCC.Core.BRepAlgoAPI import BRepAlgoAPI_Cut

import ENigMAocc
        
box = BRepPrimAPI_MakeBox(8.0, 8.0, 8.0).Shape()
cylinder = BRepPrimAPI_MakeCylinder(3.0, 6.0).Shape() 

move = gp_Trsf()
move.SetTranslation(gp_Vec(-4.0, -4.0, -4.0))
box = BRepBuilderAPI_Transform(box, move, False).Shape()

shape = BRepAlgoAPI_Cut(box, cylinder).Shape()

mesh = ENigMAocc.meshShape(shape, 0.5, 1E-3)
ENigMAocc.saveMeshFile(mesh, "occ_04a.msh")

A Fused Shape

Code

from OCC.Core.gp import gp_Vec, gp_Trsf
from OCC.Core.BRepPrimAPI import BRepPrimAPI_MakeBox, BRepPrimAPI_MakeCylinder
from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_Transform
from OCC.Core.BRepAlgoAPI import BRepAlgoAPI_Fuse

import ENigMAocc
        
box = BRepPrimAPI_MakeBox(8.0, 8.0, 8.0).Shape()
cylinder = BRepPrimAPI_MakeCylinder(3.0, 6.0).Shape() 

move = gp_Trsf()
move.SetTranslation(gp_Vec(-4.0, -4.0, -4.0))
box = BRepBuilderAPI_Transform(box, move, False).Shape()

shape = BRepAlgoAPI_Fuse(box, cylinder).Shape()

mesh = ENigMAocc.meshShape(shape, 0.5, 1E-3)
ENigMAocc.saveMeshFile(mesh, "occ_04b.msh")

More Cut-outs

Code

from OCC.Core.gp import gp_Vec, gp_Trsf
from OCC.Core.BRepPrimAPI import BRepPrimAPI_MakeBox
from OCC.Core.BRepPrimAPI import BRepPrimAPI_MakeCylinder
from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_Transform
from OCC.Core.BRepAlgoAPI import BRepAlgoAPI_Cut

import ENigMAocc
        
box = BRepPrimAPI_MakeBox(100.0, 50.0, 4.0).Shape() 
cylinder1 = BRepPrimAPI_MakeCylinder(3.5, 10.0).Shape() 

move = gp_Trsf()

cut = box

for i in range(0, 9):
    for j in range(0, 4):
        move.SetTranslation(gp_Vec(i*10 + 10, j*10 + 10, 0.0))
        cylinder_copy = BRepBuilderAPI_Transform(cylinder1, move, True).Shape()
        cut = BRepAlgoAPI_Cut(cut, cylinder_copy).Shape()
    
cylinder2 = BRepPrimAPI_MakeCylinder(6.0, 10.0).Shape()  

move.SetTranslation(gp_Vec(0, 0, 0.0))
cylinder_copy = BRepBuilderAPI_Transform(cylinder2, move, True).Shape()
cut = BRepAlgoAPI_Cut(cut, cylinder_copy).Shape()

move.SetTranslation(gp_Vec(100, 0, 0.0))
cylinder_copy = BRepBuilderAPI_Transform(cylinder2, move, True).Shape()
cut = BRepAlgoAPI_Cut(cut, cylinder_copy).Shape()

move.SetTranslation(gp_Vec(100, 50, 0.0))
cylinder_copy = BRepBuilderAPI_Transform(cylinder2, move, True).Shape()
cut = BRepAlgoAPI_Cut(cut, cylinder_copy).Shape()

move.SetTranslation(gp_Vec(0, 50, 0.0))
cylinder_copy = BRepBuilderAPI_Transform(cylinder2, move, True).Shape()
shape = BRepAlgoAPI_Cut(cut, cylinder_copy).Shape()
        
mesh = ENigMAocc.meshShape(shape, 2.0, 1E-3)
ENigMAocc.saveMeshFile(mesh, "occ_05.msh")

A STEP File

Code

from OCC.Core.STEPControl import STEPControl_Reader

import ENigMAocc

step_reader = STEPControl_Reader()
step_reader.ReadFile("lego.step")

step_reader.TransferRoot(1)
shape = step_reader.Shape(1)
    
mesh = ENigMAocc.meshShape(shape, 1.0, 1E-3)
ENigMAocc.saveMeshFile(mesh, "occ_06.msh")

The OCC Bottle

Code

import math

from OCC.Core.gp import (gp_Pnt, gp_OX, gp_Vec, gp_Trsf, gp_DZ, gp_Ax2, gp_Ax3,
                         gp_Pnt2d, gp_Dir2d, gp_Ax2d, gp_Pln)
from OCC.Core.GC import GC_MakeArcOfCircle, GC_MakeSegment
from OCC.Core.GCE2d import GCE2d_MakeSegment
from OCC.Core.Geom import Geom_CylindricalSurface
from OCC.Core.Geom2d import Geom2d_Ellipse, Geom2d_TrimmedCurve
from OCC.Core.BRepBuilderAPI import (BRepBuilderAPI_MakeEdge, BRepBuilderAPI_MakeWire,
                                     BRepBuilderAPI_MakeFace, BRepBuilderAPI_Transform)
from OCC.Core.BRepPrimAPI import BRepPrimAPI_MakePrism, BRepPrimAPI_MakeCylinder
from OCC.Core.BRepFilletAPI import BRepFilletAPI_MakeFillet
from OCC.Core.BRepAlgoAPI import BRepAlgoAPI_Fuse
from OCC.Core.BRepOffsetAPI import BRepOffsetAPI_MakeThickSolid, BRepOffsetAPI_ThruSections
from OCC.Core.BRepLib import breplib
from OCC.Core.BRep import BRep_Builder
from OCC.Core.GeomAbs import GeomAbs_Plane
from OCC.Core.BRepAdaptor import BRepAdaptor_Surface
from OCC.Core.TopoDS import topods, TopoDS_Compound, TopoDS_Face
from OCC.Core.TopExp import TopExp_Explorer
from OCC.Core.TopAbs import TopAbs_EDGE, TopAbs_FACE
from OCC.Core.TopTools import TopTools_ListOfShape

import ENigMAocc

def face_is_plane(face: TopoDS_Face) -> bool:
    """
    Returns True if the TopoDS_Face is a plane, False otherwise
    """
    surf = BRepAdaptor_Surface(face, True)
    surf_type = surf.GetType()
    return surf_type == GeomAbs_Plane


def geom_plane_from_face(aFace: TopoDS_Face) -> gp_Pln:
    """
    Returns the geometric plane entity from a planar surface
    """
    return BRepAdaptor_Surface(aFace, True).Plane()


height = 70
width = 50
thickness = 30

print("creating bottle")
# The points we'll use to create the profile of the bottle's body
aPnt1 = gp_Pnt(-width / 2.0, 0, 0)
aPnt2 = gp_Pnt(-width / 2.0, -thickness / 4.0, 0)
aPnt3 = gp_Pnt(0, -thickness / 2.0, 0)
aPnt4 = gp_Pnt(width / 2.0, -thickness / 4.0, 0)
aPnt5 = gp_Pnt(width / 2.0, 0, 0)

aArcOfCircle = GC_MakeArcOfCircle(aPnt2, aPnt3, aPnt4)
aSegment1 = GC_MakeSegment(aPnt1, aPnt2)
aSegment2 = GC_MakeSegment(aPnt4, aPnt5)

# Could also construct the line edges directly using the points instead of the resulting line
aEdge1 = BRepBuilderAPI_MakeEdge(aSegment1.Value())
aEdge2 = BRepBuilderAPI_MakeEdge(aArcOfCircle.Value())
aEdge3 = BRepBuilderAPI_MakeEdge(aSegment2.Value())

# Create a wire out of the edges
aWire = BRepBuilderAPI_MakeWire(aEdge1.Edge(), aEdge2.Edge(), aEdge3.Edge())

# Quick way to specify the X axis
xAxis = gp_OX()

# Set up the mirror
aTrsf = gp_Trsf()
aTrsf.SetMirror(xAxis)

# Apply the mirror transformation
aBRespTrsf = BRepBuilderAPI_Transform(aWire.Wire(), aTrsf)

# Get the mirrored shape back out of the transformation and convert back to a wire
aMirroredShape = aBRespTrsf.Shape()

# A wire instead of a generic shape now
aMirroredWire = topods.Wire(aMirroredShape)

# Combine the two constituent wires
mkWire = BRepBuilderAPI_MakeWire()
mkWire.Add(aWire.Wire())
mkWire.Add(aMirroredWire)
myWireProfile = mkWire.Wire()

# The face that we'll sweep to make the prism
myFaceProfile = BRepBuilderAPI_MakeFace(myWireProfile)

# We want to sweep the face along the Z axis to the height
aPrismVec = gp_Vec(0, 0, height)
myBody_step1 = BRepPrimAPI_MakePrism(myFaceProfile.Face(), aPrismVec)

# Add fillets to all edges through the explorer
mkFillet = BRepFilletAPI_MakeFillet(myBody_step1.Shape())
anEdgeExplorer = TopExp_Explorer(myBody_step1.Shape(), TopAbs_EDGE)

while anEdgeExplorer.More():
    anEdge = topods.Edge(anEdgeExplorer.Current())
    mkFillet.Add(thickness / 12.0, anEdge)

    anEdgeExplorer.Next()

# Create the neck of the bottle
neckLocation = gp_Pnt(0, 0, height)
neckAxis = gp_DZ()
neckAx2 = gp_Ax2(neckLocation, neckAxis)

myNeckRadius = thickness / 4.0
myNeckHeight = height / 10.0

mkCylinder = BRepPrimAPI_MakeCylinder(neckAx2, myNeckRadius, myNeckHeight)

myBody_step2 = BRepAlgoAPI_Fuse(mkFillet.Shape(), mkCylinder.Shape())

# Our goal is to find the highest Z face and remove it
zMax = -1.

# We have to work our way through all the faces to find the highest Z face so we can remove it for the shell
aFaceExplorer = TopExp_Explorer(myBody_step2.Shape(), TopAbs_FACE)
while aFaceExplorer.More():
    aFace = topods.Face(aFaceExplorer.Current())
    if face_is_plane(aFace):
        aPlane = geom_plane_from_face(aFace)

        # We want the highest Z face, so compare this to the previous faces
        aPntLoc = aPlane.Location()
        aZ = aPntLoc.Z()
        if aZ > zMax:
            zMax = aZ
    aFaceExplorer.Next()

facesToRemove = TopTools_ListOfShape()
facesToRemove.Append(aFace)

myBody_step3 = BRepOffsetAPI_MakeThickSolid(myBody_step2.Shape(), facesToRemove, -thickness / 50.0, 0.001)

# Set up our surfaces for the threading on the neck
neckAx2_Ax3 = gp_Ax3(neckLocation, gp_DZ())
aCyl1 = Geom_CylindricalSurface(neckAx2_Ax3, myNeckRadius * 0.99)
aCyl2 = Geom_CylindricalSurface(neckAx2_Ax3, myNeckRadius * 1.05)

# Set up the curves for the threads on the bottle's neck
aPnt = gp_Pnt2d(2.0 * math.pi, myNeckHeight / 2.0)
aDir = gp_Dir2d(2.0 * math.pi, myNeckHeight / 4.0)
anAx2d = gp_Ax2d(aPnt, aDir)

aMajor = 2.0 * math.pi
aMinor = myNeckHeight / 10.0

anEllipse1 = Geom2d_Ellipse(anAx2d, aMajor, aMinor)
anEllipse2 = Geom2d_Ellipse(anAx2d, aMajor, aMinor / 4.0)

anArc1 = Geom2d_TrimmedCurve(anEllipse1, 0, math.pi)
anArc2 = Geom2d_TrimmedCurve(anEllipse2, 0, math.pi)

anEllipsePnt1 = anEllipse1.Value(0)
anEllipsePnt2 = anEllipse1.Value(math.pi)

aSegment = GCE2d_MakeSegment(anEllipsePnt1, anEllipsePnt2)

# Build edges and wires for threading
anEdge1OnSurf1 = BRepBuilderAPI_MakeEdge(anArc1, aCyl1)
anEdge2OnSurf1 = BRepBuilderAPI_MakeEdge(aSegment.Value(), aCyl1)
anEdge1OnSurf2 = BRepBuilderAPI_MakeEdge(anArc2, aCyl2)
anEdge2OnSurf2 = BRepBuilderAPI_MakeEdge(aSegment.Value(), aCyl2)

threadingWire1 = BRepBuilderAPI_MakeWire(anEdge1OnSurf1.Edge(), anEdge2OnSurf1.Edge())
threadingWire2 = BRepBuilderAPI_MakeWire(anEdge1OnSurf2.Edge(), anEdge2OnSurf2.Edge())

# Compute the 3D representations of the edges/wires
breplib.BuildCurves3d(threadingWire1.Shape())
breplib.BuildCurves3d(threadingWire2.Shape())

# Create the surfaces of the threading
aTool = BRepOffsetAPI_ThruSections(True)
aTool.AddWire(threadingWire1.Wire())
aTool.AddWire(threadingWire2.Wire())
aTool.CheckCompatibility(False)
myThreading = aTool.Shape()

# Build the resulting compound
bottle = TopoDS_Compound()
aBuilder = BRep_Builder()
aBuilder.MakeCompound(bottle)
aBuilder.Add(bottle, myBody_step3.Shape())
aBuilder.Add(bottle, myThreading)
print("bottle finished")

mesh = ENigMAocc.meshShape(bottle, 1.5, 1E-3)
ENigMAocc.saveMeshFile(mesh, "occ_07.msh")

Structural Analysis

A Cantilever

A cantilever beam fixed at x=0 with an applied force of 1000 at x=1.

Max deflection (calculated)  = -0.021174580524850536
Max deflection (theoretical) = -0.021333333333333326

Code

import math
from ENigMA import ENigMA

b = 0.05
h = 0.05

F = -1000.0
L = 1.0
E = 30E+9
I = b * h * h * h / 12

vertex1 = ENigMA.CMshNode(0.0, 0.0, 0.0)
vertex2 = ENigMA.CMshNode(L, 0.0, 0.0)
vertex3 = ENigMA.CMshNode(L, h, 0.0)
vertex4 = ENigMA.CMshNode(0.0, h, 0.0)

quadrilateral = ENigMA.CGeoQuadrilateral()

quadrilateral.addVertex(vertex1)
quadrilateral.addVertex(vertex2)
quadrilateral.addVertex(vertex3)
quadrilateral.addVertex(vertex4)

basicMesher = ENigMA.CMshBasicMesher()

basicMesher.generate(quadrilateral, 400, 20, True)

surfaceMesh = basicMesher.mesh()

for i in range(0, surfaceMesh.nbElements()):
    elementId = surfaceMesh.elementId(i)
    surfaceMesh.element(elementId).setThickness(b)

# Material
material = ENigMA.CMatMaterial()

material.addProperty(ENigMA.PT_ELASTIC_MODULUS, E)
material.addProperty(ENigMA.PT_POISSON_COEFFICIENT, 0.3)

# Stress field
u = ENigMA.CPdeField()

u.setMesh(surfaceMesh)
u.setMaterial(material)
u.setSimulationType(ENigMA.ST_STRUCTURAL)
u.setDiscretMethod(ENigMA.DM_FEM)
u.setDiscretOrder(ENigMA.DO_LINEAR)
u.setDiscretLocation(ENigMA.DL_NODE)
u.setNbDofs(2)

index = 0

for i in range(0, surfaceMesh.nbNodes()):
    nodeId = surfaceMesh.nodeId(i)
    node = surfaceMesh.node(nodeId)
    
    if (math.fabs(node.x() - 0.0) < 1E-6):
        u.setFixedValue(surfaceMesh.nodeIndex(nodeId), 0, 0.0)
        u.setFixedValue(surfaceMesh.nodeIndex(nodeId), 1, 0.0)
        
    if (math.fabs(node.x() - L) < 1E-6 and math.fabs(node.y() - h) < 1E-6):
        index = i
        u.setSource(surfaceMesh.nodeIndex(nodeId), 1, F)

sleSystem = ENigMA.laplacian(u)
sleSystem.setRhs(0)
        
pdeEquation = ENigMA.CPdeEquation(sleSystem)

pdeEquation.setSources(u)

pdeEquation.setElimination(u)

pdeEquation.solve(u)

posGmsh = ENigMA.CPosGmsh()
posGmsh.save(u, "fem_01.msh", "Structural")

# Theoretical displacement
deflection = (F * L * L * L) / (3 * E * I)

print('Max deflection (calculated)  = ' + str(u.value(index * 2 + 1)))
print('Max deflection (theoretical) = ' + str(deflection))

Thermal Analysis

A Plate

Two-dimensional steady-state thermal analysis in a rectangular plate heated at the top (T = 1), other sides kept at 0.

Temperature at point (0.5, 0.5) (calculated)  = 0.25
Temperature at point (0.5, 0.5) (theoretical) = 0.25

Code (FEM)

import math
from ENigMA import ENigMA

vertex1 = ENigMA.CMshNode(0.0, 0.0, 0.0)
vertex2 = ENigMA.CMshNode(1.0, 0.0, 0.0)
vertex3 = ENigMA.CMshNode(1.0, 1.0, 0.0)
vertex4 = ENigMA.CMshNode(0.0, 1.0, 0.0)

quadrilateral = ENigMA.CGeoQuadrilateral()

quadrilateral.addVertex(vertex1)
quadrilateral.addVertex(vertex2)
quadrilateral.addVertex(vertex3)
quadrilateral.addVertex(vertex4)

basicMesher = ENigMA.CMshBasicMesher()

basicMesher.generate(quadrilateral, 16, 16, True)

surfaceMesh = basicMesher.mesh()

# Temperature field
u = ENigMA.CPdeField()

u.setMesh(surfaceMesh)
u.setSimulationType(ENigMA.ST_THERMAL)
u.setDiscretMethod(ENigMA.DM_FEM)
u.setDiscretOrder(ENigMA.DO_LINEAR)
u.setDiscretLocation(ENigMA.DL_NODE)
u.setNbDofs(1)

index = 0

for i in range(0, surfaceMesh.nbNodes()):
    nodeId = surfaceMesh.nodeId(i)
    node = surfaceMesh.node(nodeId)
    
    if (math.fabs(node.x() - 0.0) < 1E-6):
        u.setFixedValue(surfaceMesh.nodeIndex(nodeId), 0.0)

    if (math.fabs(node.x() - 1.0) < 1E-6):
        u.setFixedValue(surfaceMesh.nodeIndex(nodeId), 0.0)

    if (math.fabs(node.y() - 0.0) < 1E-6):
        u.setFixedValue(surfaceMesh.nodeIndex(nodeId), 0.0)
        
    if (math.fabs(node.y() - 1.0) < 1E-6):
        u.setFixedValue(surfaceMesh.nodeIndex(nodeId), 1.0)

sleSystem = ENigMA.laplacian(u)
sleSystem.setRhs(0)
        
pdeEquation = ENigMA.CPdeEquation(sleSystem)

pdeEquation.setSources(u)

pdeEquation.setElimination(u)

pdeEquation.solve(u)

posGmsh = ENigMA.CPosGmsh()
posGmsh.save(u, "fem_02.msh", "Thermal")

for i in range(0, surfaceMesh.nbNodes()):
    nodeId = surfaceMesh.nodeId(i)
    node = surfaceMesh.node(nodeId)
    
    if (math.fabs(node.x() - 0.5) < 1E-6 and math.fabs(node.y() - 0.5) < 1E-6):
        index = surfaceMesh.nodeIndex(nodeId)

anaTemperature = ENigMA.CAnaTemperature()

ut = anaTemperature.steadyStateHeatConduction2D(0.5, 0.5)

print('Temperature at point (0.5, 0.5) (theoretical) = ' + str(ut))
print('Temperature at point (0.5, 0.5) (calculated)  = ' + str(u.value(index)))

Code (FVM)

import math
from ENigMA import ENigMA

vertex1 = ENigMA.CMshNode(0.0, 0.0, 0.0)
vertex2 = ENigMA.CMshNode(1.0, 0.0, 0.0)
vertex3 = ENigMA.CMshNode(1.0, 1.0, 0.0)
vertex4 = ENigMA.CMshNode(0.0, 1.0, 0.0)
vertex5 = ENigMA.CMshNode(0.0, 0.0, 0.1)
vertex6 = ENigMA.CMshNode(1.0, 0.0, 0.1)
vertex7 = ENigMA.CMshNode(1.0, 1.0, 0.1)
vertex8 = ENigMA.CMshNode(0.0, 1.0, 0.1)

hexahedron = ENigMA.CGeoHexahedron()

hexahedron.addVertex(vertex1)
hexahedron.addVertex(vertex2)
hexahedron.addVertex(vertex3)
hexahedron.addVertex(vertex4)
hexahedron.addVertex(vertex5)
hexahedron.addVertex(vertex6)
hexahedron.addVertex(vertex7)
hexahedron.addVertex(vertex8)

basicMesher = ENigMA.CMshBasicMesher()

basicMesher.generate(hexahedron, 65, 65, 1, False)

volumeMesh = basicMesher.mesh()

volumeMesh.generateFaces(1E-6);
volumeMesh.calculateFaceCentroid();
volumeMesh.calculateElementCentroid();
    
# Temperature field
u = ENigMA.CPdeField()

u.setMesh(volumeMesh)
u.setSimulationType(ENigMA.ST_THERMAL)
u.setDiscretMethod(ENigMA.DM_FVM)
u.setDiscretOrder(ENigMA.DO_LINEAR)
u.setDiscretLocation(ENigMA.DL_ELEMENT_CENTER)
u.setNbDofs(1)

index = 0

for i in range(0, volumeMesh.nbFaces()):
    faceId = volumeMesh.faceId(i)
    face = volumeMesh.face(faceId)
    
    if (math.fabs(volumeMesh.faceCentroid(faceId).x() - 0.0) < 1E-6):
        fixedBC = ENigMA.CPdeBoundaryCondition(ENigMA.BT_GENERIC_FIXED_VALUE);
        fixedBC.addCondition(ENigMA.CT_GENERIC_FIXED_VALUE, 0.0);
        u.addBCFace(faceId, fixedBC);

    if (math.fabs(volumeMesh.faceCentroid(faceId).x() - 1.0) < 1E-6):
        fixedBC = ENigMA.CPdeBoundaryCondition(ENigMA.BT_GENERIC_FIXED_VALUE);
        fixedBC.addCondition(ENigMA.CT_GENERIC_FIXED_VALUE, 0.0);
        u.addBCFace(faceId, fixedBC);

    if (math.fabs(volumeMesh.faceCentroid(faceId).y() - 0.0) < 1E-6):
        fixedBC = ENigMA.CPdeBoundaryCondition(ENigMA.BT_GENERIC_FIXED_VALUE);
        fixedBC.addCondition(ENigMA.CT_GENERIC_FIXED_VALUE, 0.0);
        u.addBCFace(faceId, fixedBC);

    if (math.fabs(volumeMesh.faceCentroid(faceId).y() - 1.0) < 1E-6):
        fixedBC = ENigMA.CPdeBoundaryCondition(ENigMA.BT_GENERIC_FIXED_VALUE);
        fixedBC.addCondition(ENigMA.CT_GENERIC_FIXED_VALUE, 1.0);
        u.addBCFace(faceId, fixedBC);

sleSystem = ENigMA.laplacian(u)
sleSystem.setRhs(0)
        
pdeEquation = ENigMA.CPdeEquation(sleSystem)
pdeEquation.solve(u)

posGmsh = ENigMA.CPosGmsh()
posGmsh.save(u, "fvm_02.msh", "Thermal")

for i in range(0, volumeMesh.nbElements()):
    elementId = volumeMesh.elementId(i)
    
    if (math.fabs(volumeMesh.elementCentroid(elementId).x() - 0.5) < 1E-6 and math.fabs(volumeMesh.elementCentroid(elementId).y() - 0.5) < 1E-6):
        index = volumeMesh.elementIndex(elementId)

anaTemperature = ENigMA.CAnaTemperature()

ut = anaTemperature.steadyStateHeatConduction2D(0.5, 0.5)

print('Temperature at point (0.5, 0.5) (theoretical) = ' + str(ut))
print('Temperature at point (0.5, 0.5) (calculated)  = ' + str(u.value(index)))

Flow Analysis

Lid-driven cavity (Re = 1000).

Note for extracting components (v0, v1) in gmsh:

Components:
MathEval (plugin)
Expression 0: v0
Expression 1: v1

or

Magnitude:
MathEval (plugin)
Expression 0: Sqrt(v0^2+v1^2)

Code FVM

import math
from ENigMA import ENigMA

def modulus(v):
    return math.sqrt(sum(x**2 for x in v))

vertex1 = ENigMA.CGeoCoordinate(+0.00, +0.00, +0.1)
vertex2 = ENigMA.CGeoCoordinate(+1.00, +0.00, +0.1)
vertex3 = ENigMA.CGeoCoordinate(+1.00, +1.00, +0.1)
vertex4 = ENigMA.CGeoCoordinate(+0.00, +1.00, +0.1)
vertex5 = ENigMA.CGeoCoordinate(+0.00, +0.00, -0.1)
vertex6 = ENigMA.CGeoCoordinate(+1.00, +0.00, -0.1)
vertex7 = ENigMA.CGeoCoordinate(+1.00, +1.00, -0.1)
vertex8 = ENigMA.CGeoCoordinate(+0.00, +1.00, -0.1)

hexahedron = ENigMA.CGeoHexahedron()

hexahedron.addVertex(vertex1)
hexahedron.addVertex(vertex2)
hexahedron.addVertex(vertex3)
hexahedron.addVertex(vertex4)
hexahedron.addVertex(vertex5)
hexahedron.addVertex(vertex6)
hexahedron.addVertex(vertex7)
hexahedron.addVertex(vertex8)

mesher = ENigMA.CMshBasicMesher()
mesher.generate(hexahedron, 80, 80, 1, False)
mesher.mesh().generateFaces(1E-3)
mesher.mesh().calculateFaceCentroid()
mesher.mesh().calculateElementCentroid()

volumeMesh = mesher.mesh()
fvmMesh = ENigMA.CFvmMesh(volumeMesh)

U = 1.0    # lid velocity
mu = 0.001 # dynamic viscosity
rho = 1.0  # density

nu = mu / rho # kinematic viscosity 

pisoSolver = ENigMA.CFvmPisoSolver(fvmMesh)

pisoSolver.setGravity(0.0, 0.0, 0.0)
pisoSolver.setMaterialProperties(rho, mu)

faceIds = ENigMA.StdVectorInt()

faceIds.clear()
for i in range(0, fvmMesh.nbFaces()):

    faceId = fvmMesh.faceId(i)
    face = fvmMesh.face(faceId)
    face.calculateCentroid()

    if (face.centroid().x() == 0.0 or
        face.centroid().x() == 1.0 or
        face.centroid().y() == 0.0 or
        face.centroid().y() == 1.0):
        faceIds.push_back(faceId)

pisoSolver.setBoundaryVelocity(faceIds, ENigMA.BT_WALL_NO_SLIP, 0.0, 0.0, 0.0)
pisoSolver.setBoundaryPressure(faceIds, ENigMA.BT_WALL_NO_SLIP, 0.0)

faceIds.clear()
for i in range(0, fvmMesh.nbFaces()):

    faceId = fvmMesh.faceId(i)
    face = fvmMesh.face(faceId)
    face.calculateCentroid()

    if (face.centroid().y() == 1.0):
        faceIds.push_back(faceId)

pisoSolver.setBoundaryVelocity(faceIds, ENigMA.BT_WALL_NO_SLIP, U, 0.0, 0.0)

# Courant < 1
dt = U / 80		

iter = 2000

# Flow in a rectangle
for ii in range(0, iter):
    pisoSolver.iterate(dt)
    r = modulus(pisoSolver.residual())
    print('Iter = {} of {} residual = {}'.format(ii + 1, iter, r))
    if (r < 1E-6):
        break
    
u = ENigMA.CPdeField()
u.setMesh(volumeMesh)
u.setSimulationType(ENigMA.ST_FLOW)
u.setDiscretMethod(ENigMA.DM_FVM)
u.setDiscretOrder(ENigMA.DO_LINEAR)
u.setDiscretLocation(ENigMA.DL_ELEMENT_CENTER)
u.setNbDofs(2)
u.setSize(fvmMesh.nbControlVolumes() * 2)
   
for i in range(0, fvmMesh.nbControlVolumes()):

    controlVolumeId = fvmMesh.controlVolumeId(i)

    u.setValue(i * 2 + 0, pisoSolver.u(controlVolumeId))
    u.setValue(i * 2 + 1, pisoSolver.v(controlVolumeId))

posGmsh = ENigMA.CPosGmsh()
posGmsh.save(u, "fvm_lid_01.msh", "Flow")

Code FEM

import math
from ENigMA import ENigMA

edgeMesh = ENigMA.CMshMesh()

node1 = ENigMA.CMshNode(0.0, 0.0, 0.0)
node2 = ENigMA.CMshNode(1.0, 0.0, 0.0)
node3 = ENigMA.CMshNode(1.0, 1.0, 0.0)
node4 = ENigMA.CMshNode(0.0, 1.0, 0.0)

edgeMesh.addNode(1, node1)
edgeMesh.addNode(2, node2)
edgeMesh.addNode(3, node3)
edgeMesh.addNode(4, node4)

element1 = ENigMA.CMshElement(ENigMA.ET_BEAM)
element1.addNodeId(1)
element1.addNodeId(2)

element2 = ENigMA.CMshElement(ENigMA.ET_BEAM)
element2.addNodeId(2)
element2.addNodeId(3)

element3 = ENigMA.CMshElement(ENigMA.ET_BEAM)
element3.addNodeId(3)
element3.addNodeId(4)

element4 = ENigMA.CMshElement(ENigMA.ET_BEAM)
element4.addNodeId(4)
element4.addNodeId(1)

edgeMesh.addElement(1, element1)
edgeMesh.addElement(2, element2)
edgeMesh.addElement(3, element3)
edgeMesh.addElement(4, element4)

triangleMesher = ENigMA.CMshTriangleMesher()

meshSize = 0.03

edgeMesh.generateFaces(1E-3)
triangleMesher.remesh(edgeMesh, meshSize)
interiorPoints = ENigMA.StdVectorCGeoCoordinate()
triangleMesher.generate(edgeMesh, 9999, interiorPoints, meshSize, meshSize, meshSize, 1E-6)
triangleMesher.flipEdges(triangleMesher.mesh())
triangleMesher.relaxNodes(triangleMesher.mesh())

surfaceMesh = triangleMesher.mesh()

U = 1.0    # lid velocity
mu = 0.001 # dynamic viscosity
rho = 1.0  # density

nu = mu / rho # kinematic viscosit

cbsSolver = ENigMA.CFemCbsSolver2(surfaceMesh)

cbsSolver.setGravity(0.0, 0.0);
cbsSolver.setMaterialProperties(rho, mu);

index = 0

for i in range(0, surfaceMesh.nbNodes()):
    nodeId = surfaceMesh.nodeId(i)
    node = surfaceMesh.node(nodeId)
    
    if (math.fabs(node.x() - 0.0) < 1E-6):
        cbsSolver.u().setFixedValue(surfaceMesh.nodeIndex(nodeId), 0.0)
        cbsSolver.v().setFixedValue(surfaceMesh.nodeIndex(nodeId), 0.0)

    if (math.fabs(node.x() - 1.0) < 1E-6):
        cbsSolver.u().setFixedValue(surfaceMesh.nodeIndex(nodeId), 0.0)
        cbsSolver.v().setFixedValue(surfaceMesh.nodeIndex(nodeId), 0.0)

    if (math.fabs(node.y() - 0.0) < 1E-6):
        cbsSolver.u().setFixedValue(surfaceMesh.nodeIndex(nodeId), 0.0)
        cbsSolver.v().setFixedValue(surfaceMesh.nodeIndex(nodeId), 0.0)
        
    if (math.fabs(node.y() - 1.0) < 1E-6):
        cbsSolver.u().setFixedValue(surfaceMesh.nodeIndex(nodeId), 1.0)
        cbsSolver.v().setFixedValue(surfaceMesh.nodeIndex(nodeId), 0.0)

    cbsSolver.u().setValue(i, 0.0);
    cbsSolver.v().setValue(i, 0.0);

    cbsSolver.p().setValue(i, 0.0);
        
# Courant < 1
dt = 0.01

iter = 20

# Flow in a rectangle
for ii in range(0, iter):
    cbsSolver.iterate(dt)
    print('Iter = {} of {}'.format(ii + 1, iter))

u = ENigMA.CPdeField()
u.setMesh(surfaceMesh)
u.setSimulationType(ENigMA.ST_FLOW)
u.setDiscretMethod(ENigMA.DM_FEM)
u.setDiscretOrder(ENigMA.DO_LINEAR)
u.setDiscretLocation(ENigMA.DL_NODE)
u.setNbDofs(2)
u.setSize(surfaceMesh.nbNodes() * 2)
   
for i in range(0, surfaceMesh.nbNodes()):
    nodeId = surfaceMesh.nodeId(i)

    u.setValue(i * 2 + 0, cbsSolver.u().value(nodeId))
    u.setValue(i * 2 + 1, cbsSolver.v().value(nodeId))

posGmsh = ENigMA.CPosGmsh()
posGmsh.save(u, "fem_lid_01.msh", "Flow")

Name		Name	Last commit message	Last commit date
Latest commit History 597 Commits
.github/workflows		.github/workflows
images		images
trunk		trunk
.clang-format		.clang-format
.codacy.yml		.codacy.yml
.gitignore		.gitignore
.lgtm.yml		.lgtm.yml
.travis.yml		.travis.yml
BuildRelease.py		BuildRelease.py
LICENSE.txt		LICENSE.txt
README.md		README.md
build.bat		build.bat
convert.bat		convert.bat
valgrind.log		valgrind.log

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Repository files navigation

ENigMA - Extended Numerical Multiphysics Analysis

Repository

Structure

Capabilities

Dependencies

Core

Tests

Examples

Install

Usage

Quick Start

Examples

Mesh Generation

A Cylinder

A Torus

A Cut-out Sphere

A Cut-out Shape

A Fused Shape

More Cut-outs

A STEP File

The OCC Bottle

Structural Analysis

A Cantilever

Thermal Analysis

A Plate

Flow Analysis

About

Releases 20

Packages

Contributors 3

Languages

License

bjaraujo/ENigMA

Folders and files

Latest commit

History

Repository files navigation

ENigMA - Extended Numerical Multiphysics Analysis

Repository

Structure

Capabilities

Dependencies

Core

Tests

Examples

Install

Usage

Quick Start

Examples

Mesh Generation

A Cylinder

A Torus

A Cut-out Sphere

A Cut-out Shape

A Fused Shape

More Cut-outs

A STEP File

The OCC Bottle

Structural Analysis

A Cantilever

Thermal Analysis

A Plate

Flow Analysis

About

Topics

Resources

License

Stars

Watchers

Forks

Releases 20

Packages 0

Contributors 3

Languages

Packages