NYE KODEN MED BESKRIVELSER

# BESO Python Script — Endringer og kodeforklaring

## Hva ble endret fra originalskriptet?

### Endring 1: Material leses fra din modell (fmtMdb)

**Original:**
```python
mdl.Material('Material01').Elastic(((1.0, 0.3), ))
mdl.HomogeneousSolidSection('sldSec','Material01')
mdl.Material('Material02').Elastic(((0.001**3, 0.3), ))
```
Hardkoder E=1.0 og ν=0.3 — dine materialegenskaper ignoreres fullstendig.

**Endret til:**
```python
myMat = mdl.materials['Material-1']
E  = myMat.elastic.table[0][0]   # Leser E fra din modell (690000 MPa)
nu = myMat.elastic.table[0][1]   # Leser ν fra din modell (0.33)

mdl.Material('Material01').Elastic(((E, nu),))
mdl.Material('Material02').Elastic(((E * 0.001**3, nu),))
```
Nå brukes dine faktiske materialverdier (E=690 000 MPa, ν=0.33 for aluminium).

---

### Endring 2: Guard mot re-kjøring (fmtMdb)

**Lagt til:**
```python
if 'Material01' in mdl.materials.keys(): del mdl.materials['Material01']
if 'Material02' in mdl.materials.keys(): del mdl.materials['Material02']
if 'sldSec' in mdl.sections.keys(): del mdl.sections['sldSec']
if 'voidSec' in mdl.sections.keys(): del mdl.sections['voidSec']
```
Hvis du kjører scriptet to ganger på samme .cae-fil vil Abaqus krasje fordi
Material01/sldSec allerede eksisterer. Denne sjekken sletter dem først hvis de finnes.

---

### Endring 3: Python 3 fix — list() rundt dict.values() (tre steder)

Dette var den kritiske feilen som krasjet koden etter første iterasjon.

**Problem:** Abaqus 2024 bruker Python 3 der dict.values() returnerer en "view",
ikke en liste. Abaqus sin interne sum()-funksjon (importert via `from abaqus import *`)
aksepterte ikke dette.

**Feil (original):**
```python
lo, hi = min(Ae.values()), max(Ae.values())   # KRASJ
if sum(Xe.values())-tv > 0: lo = th           # KRASJ
vh.append(sum(xe.values())/len(xe))           # KRASJ
```

**Fikset:**
```python
lo, hi = min(list(Ae.values())), max(list(Ae.values()))
if sum(list(Xe.values()))-tv > 0: lo = th
vh.append(sum(list(xe.values()))/len(xe))
```
`list()` konverterer view til en vanlig Python-liste som sum()/min()/max() håndterer.

---

## Komplett kodeforklaring

### Importer (linje 1-4)
```python
import math, customKernel
from abaqus import getInput, getInputs
from odbAccess import openOdb
```
- `math` — for sqrt og fabs
- `customKernel` — Abaqus-modul for å lagre egne data i .cae-filen (brukes til å lagre konvergenshistorikk)
- `getInput/getInputs` — viser dialogbokser til brukeren
- `openOdb` — åpner Abaqus output database (.odb) for å lese resultater

---

### fmtMdb — Klargjør modellen (linje 7-27)

Denne funksjonen kjøres én gang før optimaliseringen starter.

```python
def fmtMdb(Mdb):
```

**Hva den gjør:**
1. Leser E og ν fra Material-1
2. Lager to materialer:
   - `Material01` (solid): dine faktiske materialverdier
   - `Material02` (void): E_void = E × 0.001³ ≈ 10⁻⁹ × E (nesten ingenting)
3. Lager to seksjoner: `sldSec` og `voidSec`
4. Starter med å tildele alle elementer til `sldSec` (fullt design)
5. Ber Abaqus om å lagre to output-variabler:
   - `ELEDEN` — element strain energy density (brukes til sensitivitet)
   - `ALLWK` — total ekstern arbeid = compliance (objektfunksjon)

**Hvorfor soft-kill?**
Void-elementer fjernes ikke fysisk, de får bare ekstremt lav stivhet.
Dette gir en glattere optimalisering og unngår singulære stivhetsmatriser.
E_void = E × xmin^p = E × 0.001³ (SIMP-modell med p=3)

---

### FEA — Kjør analyse og hent sensitiviteter (linje 30-38)

```python
def FEA(Iter, Mdb, Xe, Ae):
```

**Hva den gjør per iterasjon:**
1. Lager og kjører en ny Abaqus-job: `Design_Job0`, `Design_Job1`, ...
2. Venter til jobben er ferdig
3. Åpner .odb-filen
4. Leser `ESEDEN` (element strain energy density) for hvert element
5. Beregner sensitivitet: α_e = E_e / x_e
6. Leser total compliance fra `ALLWK`
7. Returnerer compliance som objektfunksjonsverdi

**Sensitiviteten α_e** forteller hvor mye compliance øker hvis elementet fjernes.
Høy α_e = elementet er viktig → behold.
Lav α_e = elementet bidrar lite → kandidat for fjerning.

---

### preFlt — Forbered filterkart (linje 41-56)

```python
def preFlt(Rmin, Elmts, Nds, Fm):
```

Kjøres kun én gang. Beregner for hvert element hvilke naboelementer som
ligger innenfor filterradius Rmin, og vektfaktorene gj for disse.

**Vektfunksjon:** w(r_ej) = max(0, Rmin - r_ej)
Elementer nærmere sentrum vektes høyere.

Resultatet lagres i en dictionary `Fm`:
```
Fm = {
  element_label: [[nabo_label1, nabo_label2, ...], [vekt1, vekt2, ...]],
  ...
}
```

**Hvorfor filtrering?**
Uten filter gir BESO sjakkmønsterinstabilitet — annenhver element
solid/void i et rutemønster som ser optimalt ut matematisk men ikke fysisk.
Filteret glatter sensitivitetene slik at nærliggende elementer påvirker hverandre.

---

### fltAe — Anvend filter (linje 59-63)

```python
def fltAe(Ae, Fm):
```

Bruker filterkartet fra preFlt til å glatte sensitivitetene:

α̃_e = Σ_j (g_j × α_j)

Der summen går over alle naboelementer j innenfor Rmin.
Hvert elements sensitivitet erstattes av et vektet gjennomsnitt av naboene.

---

### BESO — Oppdater designvariabler (linje 66-81)

```python
def BESO(Vf, Xe, Ae, Part, Elmts):
```

**Hva den gjør:**
1. Finner terskelverdien th via biseksjon slik at riktig antall elementer
   blir solid og vi treffer målvolumet Vf
2. Setter Xe[e] = 1.0 (solid) hvis α_e > th, ellers 0.001 (void)
3. Deler alle elementer i to lister: slb (solid) og vlb (void)
4. Tildeler seksjonene sldSec/voidSec til de to element-settene i Abaqus

**Biseksjonsalgoritmen:**
```
lo = min(sensitivitet), hi = max(sensitivitet)
Gjenta:
    th = (lo + hi) / 2
    Tell elementer med α_e > th
    Hvis for mange solide: lo = th
    Hvis for få solide:    hi = th
Inntil (hi-lo)/hi < 0.00001
```

---

### Hovedprogram (linje 84-115)

**Blokk 1: Input**
```python
pars = (('VolFrac:', '0.5'), ('Rmin:', '1'), ('ER:', '0.02'))
vf, rmin, ert = [float(k) ... for k in getInputs(pars)]
mddb = openMdb(getInput('Input CAE file:', default='Testing.cae'))
```
Henter parametere og åpner .cae-filen din.

**Blokk 2: Initialisering**
```python
fmtMdb(mddb)
for el in elmts: xe[el.label] = 1.0   # Start med fullt design
if rmin > 0: preFlt(rmin, elmts, nds, fm)
```
Alle elementer starter som solide (xe=1.0). Filter forberedes.

**Blokk 3: Iterasjonsloop**
```python
while change > 0.001:
    iter += 1
    oh.append(FEA(iter, mddb, xe, ae))          # 1. Kjør FEA
    if rmin > 0: fltAe(ae, fm)                  # 2. Filtrer sensitiviteter
    if iter > 0: ae = dict([(k,(ae[k]+oae[k])/2.0) ...])  # 3. Historisk gjennomsnitt
    oae = ae.copy()
    vh.append(sum(list(xe.values()))/len(xe))   # 4. Registrer volum
    nv = max(vf, vh[-1]*(1.0-ert))              # 5. Beregn målvolum
    BESO(nv, xe, ae, part, elmts)               # 6. Oppdater topologi
    if iter > 10: change = ...                  # 7. Sjekk konvergens
```

Hvert steg i loopen:
1. FEA med nåværende topologi → compliance og sensitiviteter
2. Glatt sensitivitetene med filterkartet
3. Gjennomsnittliggjør med forrige iterasjon (stabiliserer konvergens)
4. Registrer nåværende volumfraksjon
5. Beregn neste målvolum: nv = max(vf, V_current × (1 - ER))
   → Volumet synker med ER=2% per iterasjon inntil det når VolFrac
6. BESO oppdaterer hvilke elementer som er solid/void
7. Etter 10 iterasjoner: sjekk om compliance har stabilisert seg

**Konvergenskriterium:**
```python
change = |Σ(C[k-4:k+1]) - Σ(C[k-9:k-4])| / Σ(C[k-9:k-4])
```
Sammenligner gjennomsnittlig compliance for siste 5 vs forrige 5 iterasjoner.
Stopper når endringen er under 0.1%.

**Blokk 4: Lagre**
```python
mddb.customData.History = {'vol': vh, 'obj': oh}
mddb.saveAs('Final_design.cae')
```
Lagrer konvergenshistorikken og den ferdige topologien.

---

## Oppsummert flyt

```
Testing.cae (ditt originale design)
        ↓
fmtMdb: Les Material-1, lag solid/void seksjoner
        ↓
preFlt: Beregn filterkart (én gang)
        ↓
[For hver iterasjon:]
    FEA → sensitiviteter og compliance
    Filtrer + historisk gjennomsnitt
    BESO: fjern elementer med lav sensitivitet
    Sjekk konvergens
        ↓
Final_design.cae (optimalisert topologi)
```