Skip to content
Snippets Groups Projects
Commit bb5694ba authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

submarine moved to ../3d_large

parent a707dead
No related branches found
No related tags found
No related merge requests found
Merge "Submarine.stp";
//Mesh.CharacteristicLengthExtendFromBoundary=0;
Mesh.Algorithm=9;
Mesh.Algorithm3D=9;
Mesh.Smoothing=0;
Mesh.Recombine3DAll=1;
//0: hex, 1: hex+prisms, 2: hex+prism+1-step pyramids, 3: hex+prism+2-steps pyramids
Mesh.Recombine3DLevel = 2;
//conformity - 0: nonconforming, 1: trihedra, 2: pyramids+trihedra, 3:pyramids+hexPrismSplit+trihedra, 4:hexPrismSplit+trihedra
Mesh.Recombine3DConformity = 4;
Mesh.SaveParametric = 1;
MeshAlgorithm Surface {12} = 5; //Not in final mesh
MeshAlgorithm Surface {10} = 5;
MeshAlgorithm Surface {11} = 5;
//Mesh.Algorithm=5;
//Mesh.Algorithm3D=5;
Mesh.Optimize = 0;
//Mesh.CharacteristicLengthMax = 25;
LcMin = 39;
LcMax = 78;
//LcMin = 80;
//LcMax = 160;
Characteristic Length {14, 13, 15, 16, 10, 9, 11, 12, 6, 5, 19, 8, 4, 3, 7, 2, 1} = LcMax;
Characteristic Length {21} = LcMin;
Characteristic Length {17, 18} = LcMin;
Characteristic Length {20} = LcMax;
Compound Line(34) = {22, 24, 25};
Compound Surface(13) = {10};
Compound Surface(14) = {11};
Compound Surface(15) = {12};
submarineLine[] = {15, 21, 17, 26, 11, 20, 13, 27, 29, 8, 19, 10, 34};
Physical Line(1) = submarineLine[];
submarineSurf[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 13};
Physical Surface(2) = submarineSurf[];
Physical Surface(3) = {14,15};
Physical Volume(4) = {1};
Physical Line(5) = {30, 31};
Source diff could not be displayed: it is too large. Options to address this: view the blob.
from gmshpy import *
import sys
import pickle
import preprofix
name = "Submarine"
g = GModel()
g.load(name + ".geo")
g.mesh(3)
g.writeMSH(name + ".msh", 3.0)
#g.load(name + ".msh")
g.setAllVolumesPositiveTopology()
g.writeMSH(name + "_topology.msh", 3.0)
preprofix.fixTetras(g)
preprofix.fixPyramids(g)
preprofix.fixPrisms(g)
preprofix.fixHexs(g)
g.writeMSH(name + "_FIX.msh", 3.0)
#exit()
#preprofix.fixPyramids(g)
#g.writeMSH(name + "Fix.msh",3.0)
OH = MeshQualOptParameters()
OH.onlyVisible = False
OH.dim = 3
OH.fixBndNodes = False # Fix boundary nodes or not
OH.strategy = 0 # 0 = Connected blobs, 1 = Adaptive one-by-one (recommended in 3D)
OH.excludeHex = False
OH.excludePrism = False
OH.nbLayers = 3 # Nb. of layers around invalid element to construct blob
OH.distanceFactor = 3 # Distance factor to construct blob
OH.maxPatchAdapt = 5 # Number of blob adaption iteration
OH.maxLayersAdaptFact = 3 # Factor to multiply nb. of layers when adapting
OH.distanceAdaptFact = 3 # Factor to multiply distance factor when adapting
#OH.weight = 0
OH.onlyValidity = False
#OH.minTargetIdealJac = 0.1
OH.minTargetInvCondNum = 0.3
#OH.weightFixed = 1.e-3
#OH.weightFree = 1.e-6
OH.nCurses = 1
OH.logFileName = "log"
OH.maxOptIter = 20 # Nb of optimixation iterations
OH.maxBarrierUpdates = 15 # Nb. of optimization passes
#print("minTargetIdealJac = %g" % OH.minTargetIdealJac)
print("minTargetInvCondNum = %g" % OH.minTargetInvCondNum)
print("nbLayers = %g, distanceFactor = %g" % (OH.nbLayers, OH.distanceFactor))
print("itMax = %g, optPassMax = %g" % (OH.maxOptIter, OH.maxBarrierUpdates))
#print("weightFree = %g, weightFixed = %g" % (OH.weightFree, OH.weightFixed))
print("strategy = %g" % OH.strategy)
print("maxAdaptBlob = %g" % OH.maxPatchAdapt)
print("maxLayersAdaptFact = %g, distanceAdaptFact = %g" % (OH.maxLayersAdaptFact, OH.distanceAdaptFact))
MeshQualityOptimizer(g, OH)
#print("RESULT: minNCJ = %g, maxNCJ = %g, CPU = %g" % (OH.minNCJ, OH.maxNCJ, OH.CPU))
#print("RESULT: minInvIdealJac = %g, maxInvIdealJac = %g, CPU = %g" %
# (OH.minIdealJac, OH.maxIdealJac, OH.CPU))
print("RESULT: minInvCondNum = %g, maxInvCondNum = %g, CPU = %g" %
(OH.minInvCondNum, OH.maxInvCondNum, OH.CPU))
g.writeMSH(name + "_opt.msh", 3.0)
pOpt = meshPartitionOptions()
pOpt.setNumOfPartitions(5)
PartitionMesh(g, pOpt)
g.writeMSH(name + "_opt_part.msh", 3.0)
#g.writeMSH(name + "_OK.msh",3.0,False,False,False,1.0,0,0)
import gmshpy
import math
def areCoplanarFaces(f0, f1, tol) :
n0 = f0.normal()
n1 = f1.normal()
return (abs(math.acos(gmshpy.dot(n0, n1))*180/math.pi) < tol)
def nbVerticesOnSurface(el) :
vent = [el.getVertex(i).onWhat().dim() for i in range(el.getNumVertices())]
return len([e for e in vent if e < 3])
def isBadElement(el) :
fs = [el.getFace(i) for i in range(el.getNumFaces())]
bad = False
for i in range(len(fs)) :
for j in range(i):
bad |= areCoplanarFaces(fs[i], fs[j], 10) #10 degrees min angle allowed
return bad
def fixTetras(m):
for r in m.bindingsGetRegions() :
goodt = [t for t in r.tetrahedra if (not ((nbVerticesOnSurface(t) == 4) and isBadElement(t)))]
oldSize = r.tetrahedra.size()
r.tetrahedra.clear()
for t in goodt :
r.addTetrahedron(t)
print("-- Removed %i tetrahedra over %i"%(oldSize-r.tetrahedra.size(), oldSize))
def fixPyramids(m):
for r in m.bindingsGetRegions() :
goodp = [p for p in r.pyramids if (not ((nbVerticesOnSurface(p) == 5) and isBadElement(p)))]
oldSize = r.pyramids.size()
r.pyramids.clear()
for p in goodp :
r.addPyramid(p)
print("-- Removed %i pyramids over %i"%(oldSize-r.pyramids.size(), oldSize))
def fixPrisms(m):
for r in m.bindingsGetRegions() :
goodp = [p for p in r.prisms if (not ((nbVerticesOnSurface(p) == 6) and isBadElement(p)))]
oldSize = r.prisms.size()
r.prisms.clear()
for p in goodp :
r.addPrism(p)
print("-- Removed %i prisms over %i"%(oldSize-r.prisms.size(), oldSize))
def fixHexs(m):
for r in m.bindingsGetRegions() :
goodh = [h for h in r.hexahedra if (not ((nbVerticesOnSurface(h) == 8) and isBadElement(p)))]
oldSize = r.hexahedra.size()
r.hexahedra.clear()
for h in goodh :
r.addHexahedron(h)
print("-- Removed %i hexahedra over %i"%(oldSize-r.hexahedra.size(), oldSize))
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment