Usage

Installation

STLCutters is a registered package. You can install it by running:

# Use ] to enter the Pkg REPL mode
pkg> add STLCutters

Load package

Load the package normally with

using STLCutters

Downloading STL files

The STLCutters package works with STL files. Therefore, first, you need to create or download an STL file. The Thingi10k dataset is one of the options for testing. You can download any indexed geometry with download_thingi10k. E.g.,

filename = download_thingi10k(293137)

Discretization

Since STLCutters is an extension of GridapEmbedded, it utilizes the same workflow to solve PDEs on embedded domains. In particular, STLCutters extends the cut function from GridapEmbedded.

We load the STL file with an STLGeometry object. E.g.,

filename = download_thingi10k(293137)
geo = STLGeometry(filename)

Then, we define a DiscreteModel around, e.g., a CartesianDiscreteModel.

pmin,pmax = get_bounding_box(geo)
model = CartesianDiscreteModel(pmin,pmax,())

Now, we can cut the model with the STLGeometry.

cutgeo = cut(model,geo)
Note

Some STL files may not be properly defined. It is recommended to check the requisites before proceeding with check_requisites. In particular, one checks if the STL is a closed surface, a manifold and has a bounded density of faces.

Usage with Gridap

Once the geometry is discretized, one can generate the embedded triangulations to solve PDEs with Gridap, see also Gridap Tutorials.

Like in GridapEmbedded, we extract the embedded triangulations as follows.

Ωact = Triangulation(cutgeo,ACTIVE)
Ω = Triangulation(cutgeo)
Γ = EmbeddedBoundary(cutgeo)
Λ = SkeletonTriangulation(cutgeo)

Serial example

Now, we provide an example of the solution of a Poisson problem on the embedded domain.

using STLCutters
using Gridap
using GridapEmbedded
cells = (10,10,10)
filename = "stl_file_path.stl"
# Domain and discretization
geo = STLGeometry(filename)
pmin,pmax = get_bounding_box(geo)
model = CartesianDiscreteModel(pmin,pmax,cells)
cutgeo = cut(model,geo)
# Cell aggregation
aggregates = aggregate(AggregateAllCutCells(),cutgeo)
# Triangulations
Ω_act = Triangulation(cutgeo,ACTIVE)
Ω = Triangulation(cutgeo)
Γ = EmbeddedBoundary(cutgeo)
nΓ = get_normal_vector(Γ)   
dΩ = Measure(Ω,2)
dΓ = Measure(Γ,2)
# FE spaces
Vstd = TestFESpace(Ω_act,ReferenceFE(lagrangian,Float64,1))
V = AgFEMSpace(Vstd,aggregates)
U = TrialFESpace(V)
# Weak form
γ = 10.0
h = (pmax - pmin)[1] / cells[1]
ud(x) = x[1] - x[2]
f = 0
a(u,v) =
    ∫( ∇(v)⋅∇(u) )dΩ +
    ∫( (γ/h)*v*u  - v*(nΓ⋅∇(u)) - (nΓ⋅∇(v))*u )dΓ
l(v) =
    ∫( v*f )dΩ +
    ∫( (γ/h)*v*ud - (nΓ⋅∇(v))*ud )dΓ
# Solve
op = AffineFEOperator(a,l,U,V)
uh = solve(op)
writevtk(Ω,"results",cellfields=["uh"=>uh])
Note

The STL file can be downloaded using download_thingi10k.

Note

One can consider a different stabilization of the small cut-cell problem instead of AgFEM. Then, the aggregate and AgFEMSpace need to be removed. See more examples in GridapEmbedded