Commit 14996f01 authored by Amira Abdel-Rahman's avatar Amira Abdel-Rahman
Browse files

multiscale madcat

parent e02a5767
Pipeline #12101 passed with stage
in 21 seconds
......@@ -10,7 +10,8 @@ CAD
tutorial
live
.DS_Store
.gif
.png
.mp4
img
\ No newline at end of file
*.gif
*.png
*.mp4
img
res3d
\ No newline at end of file
......@@ -32,7 +32,7 @@
var interaction= {
exaggeration:1,
speed:100
speed:100*2 //100
}
var gui = new dat.GUI();
......@@ -125,6 +125,7 @@
new THREE.MeshLambertMaterial({
color: getNodeColor(restrained,loaded,displaced),
transparent: opacity>0.0?false:true,
transparent:false,
opacity: opacity
})
))
......@@ -340,12 +341,14 @@
new THREE.MeshLambertMaterial({
color: getNodeColor(restrained,loaded,displaced),
transparent: opacity>0.0?false:true,
transparent: false,
opacity: opacity
})
);
m.rotation.x=ax;
m.rotation.y=ay;
m.rotation.z=az;
m.visible = !(opacity>0.0?false:true);
return m;
......
This diff is collapsed.
This diff is collapsed.
# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2020
######################### 1. Voxel Design ###########################
setup = Dict()
### 1.b Draw Lattice
voxelSize=0.001
latticeSizeX=1
latticeSizeY=1
latticeSizeZ=1
gridSize=50
setup["latticeSizeX"]=latticeSizeX
setup["latticeSizeY"]=latticeSizeY
setup["latticeSizeZ"]=latticeSizeZ
setup["gridSize"]=gridSize
setup["rhino"]=false;
setup["useVoxelList"]=true;
######################### 2. Boundary Conditions #########################
######################### 2.a. Global Settings #########################
#scaling params
setup["voxelSize"]=voxelSize; #voxel size
setup["scale"]=1e4; #scale for visualization
setup["hierarchical"]=true; #hierachical simualtion
setup["multiscale"]=true; #multiscale simualtion
#simulation params
setup["numTimeSteps"]=1000; #num of saved timesteps for simulation
setup["poisson"]=false; # account for poisson ration (only for hierarchical)
setup["linear"]=true; # linear vs non-linear
setup["thermal"]=false; #if there is change in temperature
setup["globalDamping"]=0.2; # (usually from 0.1 to 0.4)
#visualization params
setup["maxNumFiles"]=10; #num of saved timesteps for visualization
######################### 2.b. Materials #########################
#default material
material1= Dict()
material1["area"]=voxelSize*voxelSize
material1["density"]=1e3
material1["stiffness"]=1e6
material1["poissonRatio"]=0.0
material1["cTE"]=0.0 #coefficient of thermal expansion
material1["scale"]=1.0 #for multiscale simulation
material2= Dict()
material2["area"]=voxelSize*voxelSize
material2["density"]=1e3
material2["stiffness"]=1e6*2.0
material2["poissonRatio"]=0.0
material2["cTE"]=0.0 #coefficient of thermal expansion
material2["scale"]=2.0 #for multiscale simulation
#//allowed x, yand z are from -gridSize to +gridSize (floor is at 0)
setup["voxelList"]=[];
for vox in voxs
if vox[2][1]==1
append!(setup["voxelList"], [[[vox[1][1],vox[1][3],vox[1][2]],material1]])
elseif vox[2][1]==2
append!(setup["voxelList"], [[[vox[1][1],vox[1][3],vox[1][2]],material2]])
end
end
###############################################################
setup["materials"]=[
# [boundingBoxMaterial1,material1],
# [boundingBoxMaterial2,material2]
];
######################### 2.c. Supports #########################
#x,y,z,rx,ry,rz (default is pinned joing i.e [false,false,false,true,true,true])
dof=[true,true,true,true,true,true]
boundingBoxSupport1=Dict()
boundingBoxSupport1["min"]=Dict()
boundingBoxSupport1["max"]=Dict()
boundingBoxSupport1["min"]["x"]= 0;
boundingBoxSupport1["min"]["y"]= 0;
boundingBoxSupport1["min"]["z"]= 0;
boundingBoxSupport1["max"]["x"]= voxelSize*2;
boundingBoxSupport1["max"]["y"]= voxelSize*10;
boundingBoxSupport1["max"]["z"]= voxelSize*gridSize;
setup["supports"]=[
[boundingBoxSupport1,dof]
];
######################### 2.d. Loads #########################
#### 2.d.1 Static Loads
load1=Dict()
load1["x"]=0.0
load1["y"]=-0.0215
load1["z"]=0.0
boundingBoxLoad1=Dict()
boundingBoxLoad1["min"]=Dict()
boundingBoxLoad1["max"]=Dict()
boundingBoxLoad1["min"]["x"]=0;
boundingBoxLoad1["min"]["y"]=0;
boundingBoxLoad1["min"]["z"]=0;
boundingBoxLoad1["max"]["x"]=voxelSize*gridSize;
boundingBoxLoad1["max"]["y"]=voxelSize*gridSize;
boundingBoxLoad1["max"]["z"]=voxelSize*gridSize;
setup["loads"]=[
[boundingBoxLoad1,load1],
];
setup["fixedDisplacements"]=[
];
#### 2.d.2 Dynamic Loads
function floorEnabled()
return false
end
function gravityEnabled()
return false
end
function externalDisplacement(currentTimeStep,N_position,N_fixedDisplacement)
return N_fixedDisplacement
end
function externalForce(currentTimeStep,N_position,N_currentPosition,N_force)
return N_force
end
function updateTemperature(currentRestLength,currentTimeStep,mat)
return currentRestLength
end
# JULIA 1.2.0
using EnhancedGJK ,GeometryTypes,LinearAlgebra,BenchmarkTools;
# make sure MeshIO is old version (0.3.2) to use geometry types and not gemoetry basics
using MeshIO ,FileIO,CoordinateTransformations,MeshCat,RegionTrees;
# using GeometryBasics
import StaticArrays: SVector
import JSON
include("../../include/asdf/adaptive_distance_fields.jl")
using .AdaptivelySampledDistanceFields
include("../../include/asdf/reference_distance.jl")
include("../../include/asdf/asdf_functions.jl")
######################### Import Mesh and SDF Functions ###########################
wing =load("../julia/examples/CAD_Rhino/wing_rot.stl")
## if convex mesh
# origin = SVector(-10., -3., -12)
# widths = SVector(24., 24, 24)
# insideOnly=true
# maxDivision=5 #how many divisions
# minDivision=6 #how many divisions
# minSize=widths[1]/(2^minDivision) #min voxel size
# maxSize=widths[1]/(2^maxDivision) #max voxel size
# atol=minSize*widths[1]/2*0.1
# println("Min Vox Size=$(minSize),Max Vox Size=$(maxSize), atol=$(atol)")
# rtol=0.0
# adaptive = AdaptivelySampledDistanceFields.ASDF(sWing, origin, widths, tol1, tol2)
# boxmesh=getVisQuadtree(adaptive,true, sdfWing)
# save("../CAD/WingMultiscale.stl", boxmesh)
## if concave mesh
resolution=128*2
mesh=wing
empty=false
voxels,verts_min,verts_max=_voxelize(mesh, resolution,empty);
###########################
origin = SVector(-10., -4., -12)
widths = SVector(24., 24, 24) #should all be the same
w=widths[1]
insideOnly=true
maxDivision=5 #how many divisions (biggest)
minDivision=6 #how many divisions (smallest)
minSize=w/(2^minDivision) #min voxel size
maxSize=w/(2^maxDivision) #max voxel size
atol=minSize*w/2*0.1
rtol=0.0
println("Desired Min Vox Size=$(minSize),Max Vox Size=$(maxSize), atol=$(atol)")
#######
# s_wing=ReferenceDistance.signed_distance(wing)
sdfWing=signed_distance1(wing)
sWing=adjustedSignedDistance(wing) #use this when convex mesh
sWingConcave=signedDistanceConcave(wing,voxels,resolution,verts_min,verts_max)
sExWingConcave=exteriorDistanceConcave(wing,voxels,resolution,verts_min,verts_max)
#####
adaptive = AdaptivelySampledDistanceFields.ASDF(sExWingConcave, origin, widths, rtol, atol,maxSize)
# adaptive = AdaptivelySampledDistanceFields.ASDF(sWing, origin, widths, rtol, atol,maxSize)
# adaptive = AdaptivelySampledDistanceFields.ASDF(sWingConcave, origin, widths, rtol, atol,maxSize)
# boxmesh=getVisQuadtreeNonConvex(adaptive,insideOnly,voxels,resolution,verts_min,verts_max)
# save("../CAD/WingMultiscale1.stl", boxmesh)
######################### Export Voxels/cubes ###########################
cubes=[]
for leaf in allleaves(adaptive.root)
for face in RegionTrees.faces(leaf.boundary)
p=(RegionTrees.center(leaf.boundary))
pointX=Int(round(p[1],digits=0)); pointY=Int(round(p[2],digits=0)); pointZ=Int(round(p[3],digits=0))
if pointInsideVoxelGrid([p[1] p[2] p[3]],voxels,resolution,verts_min,verts_max) #println(leaf.boundary.widths)
append!(cubes,[leaf.boundary])
end
end
end
cubes=unique(cubes)
sort!(cubes, by = x -> x.widths[2]);
minSize=cubes[1].widths[1]
maxSize=cubes[end].widths[1]
println("Real Min Vox Size=$(minSize),Max Vox Size=$(maxSize)")
sort!(cubes, by = x -> x.origin[1]);
minOrigin1=cubes[1].origin[1]
sort!(cubes, by = x -> x.origin[2]);
minOrigin2=cubes[1].origin[2]
sort!(cubes, by = x -> x.origin[3]);
minOrigin3=cubes[1].origin[3]
minOrigin=[minOrigin1,minOrigin2,minOrigin3].+minSize/2.0
voxs=[]
for cube in cubes
size=cube.widths[1]/minSize
append!(voxs,[[((cube.origin.-minOrigin)./minSize .+(0.5*size)) ,[size]]])
end
module AdaptivelySampledDistanceFields
using StaticArrays
using RegionTrees
import RegionTrees: needs_refinement, refine_data
using Interpolations
@generated function evaluate(itp::AbstractInterpolation, point::SVector{N}) where N
Expr(:call, :itp, [:(point[$i]) for i in 1:N]...)
end
function evaluate(itp::AbstractInterpolation, point::AbstractArray)
itp(point...)
end
function evaluate(cell::Cell{D}, point::AbstractArray) where D <: AbstractInterpolation
leaf = findleaf(cell, point)
evaluate(leaf.data, leaf.boundary, point)
end
function evaluate(interp::AbstractInterpolation, boundary::HyperRectangle, point::AbstractArray)
coords = (point .- boundary.origin) ./ boundary.widths .+ 1
evaluate(interp, coords)
end
mutable struct SignedDistanceRefinery{F <: Function} <: AbstractRefinery
signed_distance_func::F
atol::Float64
rtol::Float64
maxSize::Float64
end
# function needs_refinement(refinery::SignedDistanceRefinery, cell::Cell)
# (maximum(cell.boundary.widths) > refinery.maxSize)||(minimum(cell.boundary.widths) > refinery.atol && needs_refinement(cell, refinery.signed_distance_func, refinery.atol, refinery.rtol))
# end
function needs_refinement(refinery::SignedDistanceRefinery, cell::Cell)
if maximum(cell.boundary.widths) > refinery.maxSize
return true
elseif minimum(cell.boundary.widths) < refinery.atol
return false
else
return (minimum(cell.boundary.widths) > refinery.atol && needs_refinement(cell, refinery.signed_distance_func, refinery.atol, refinery.rtol))
end
end
function needs_refinement(cell::Cell, signed_distance_func, atol, rtol)
for c in body_and_face_centers(cell.boundary)
value_interp = evaluate(cell, c)
value_true = signed_distance_func(c)
if !isapprox(value_interp, value_true, rtol=rtol, atol=atol)
return true
end
end
false
end
function refine_data(refinery::SignedDistanceRefinery, cell::Cell, indices)
refine_data(refinery, child_boundary(cell, indices))
end
function refine_data(refinery::SignedDistanceRefinery, boundary::HyperRectangle)
interpolate!(refinery.signed_distance_func.(vertices(boundary)),
BSpline(Linear()))
end
struct ASDF{C <: Cell} <: Function
root::C
end
function ASDF(signed_distance::Function, origin::AbstractArray,
widths::AbstractArray,
rtol=1e-2,
atol=1e-2,
maxSize=Inf64)
refinery = SignedDistanceRefinery(signed_distance, atol, rtol,maxSize)
boundary = HyperRectangle(origin, widths)
root = Cell(boundary, refine_data(refinery, boundary))
adaptivesampling!(root, refinery)
ASDF(root)
end
function (field::ASDF)(point::AbstractArray)
evaluate(field.root, point)
end
end
function getVisQuadtree(adaptive,inside,trueSDF)
verts = GeometryTypes.Point{3, Float64}[]
faces = GeometryTypes.Face{3, Int}[]
count=0
countInside=0
for leaf in allleaves(adaptive.root)
for face in RegionTrees.faces(leaf.boundary)
if (inside)
if (trueSDF(RegionTrees.center(leaf.boundary))<0.0)
i = length(verts)
for vert in face
push!(verts, GeometryTypes.Point(vert[1], vert[2], vert[3]))
end
push!(faces, GeometryTypes.Face(i+1, i+2, i+4))
push!(faces, GeometryTypes.Face(i+4, i+3, i+1))
countInside+=1
end
else
i = length(verts)
for vert in face
push!(verts, GeometryTypes.Point(vert[1], vert[2], vert[3]))
end
push!(faces, GeometryTypes.Face(i+1, i+2, i+4))
push!(faces, GeometryTypes.Face(i+4, i+3, i+1))
count+=1
end
end
end
if inside
print("inside:")
println(countInside)
else
print("all:")
println(count)
end
return GeometryTypes.HomogenousMesh(verts, faces)
end
function getVisQuadtreeNonConvex(adaptive,inside,voxels,resolution,verts_min,verts_max)
verts = GeometryTypes.Point{3, Float64}[]
faces = GeometryTypes.Face{3, Int}[]
count=0
countInside=0
for leaf in allleaves(adaptive.root)
for face in RegionTrees.faces(leaf.boundary)
if (inside)
p=(RegionTrees.center(leaf.boundary))
pointX=Int(round(p[1],digits=0)); pointY=Int(round(p[2],digits=0)); pointZ=Int(round(p[3],digits=0))
if pointInsideVoxelGrid([p[1] p[2] p[3]],voxels,resolution,verts_min,verts_max)
# if pointInsideVoxelGrid([pointX pointY pointZ],voxels,resolution,verts_min,verts_max)
i = length(verts)
for vert in face
push!(verts, GeometryTypes.Point(vert[1], vert[2], vert[3]))
end
push!(faces, GeometryTypes.Face(i+1, i+2, i+4))
push!(faces, GeometryTypes.Face(i+4, i+3, i+1))
countInside+=1
end
else
i = length(verts)
for vert in face
push!(verts, GeometryTypes.Point(vert[1], vert[2], vert[3]))
end
push!(faces, GeometryTypes.Face(i+1, i+2, i+4))
push!(faces, GeometryTypes.Face(i+4, i+3, i+1))
count+=1
end
end
end
if inside
print("inside:")
println(countInside)
else
print("all:")
println(count)
end
return GeometryTypes.HomogenousMesh(verts, faces)
end
function getHyperCubes(adaptive,trueSDF)
cubes=[]
count=0
for leaf in allleaves(adaptive.root)
for face in RegionTrees.faces(leaf.boundary)
if(trueSDF(RegionTrees.center(leaf.boundary))<0.0)# if inside
#println(leaf.boundary.widths)
append!(cubes,[leaf.boundary])
end
count+=1
end
end
cubes=unique(cubes)
println(count)
return cubes
end
function getOrderedBins(adaptive,trueSDF)
function getSDF(x)
trueSDF(RegionTrees.center(x))
end
cubes=[]
count=0
for leaf in allleaves(adaptive.root)
for face in RegionTrees.faces(leaf.boundary)
if(trueSDF(RegionTrees.center(leaf.boundary))<0.0)# if inside
#println(leaf.boundary.widths)
append!(cubes,[leaf.boundary])
end
count+=1
end
end
cubes=unique(cubes)
sort!(cubes, by = x -> x.origin[2]);
#put in bins by y value (down first)
orderedY=[]
currentY=cubes[1].origin[2];
currentLayer=1
append!(orderedY,[[]]);
for cube in cubes
if(cube.origin[2]==currentY)
append!(orderedY[currentLayer],[cube]);
else
append!(orderedY,[[]]);
currentY=cube.origin[2]
currentLayer=currentLayer+1
append!(orderedY[currentLayer],[cube]);
end
end
#for each orderY bin put in bins according to size
orderedYandSize=[]
count=1
for binY in orderedY
append!(orderedYandSize,[[]]);
orderedY[count]=sort(binY, rev = true,by = x -> x.widths[1])
count+=1
end
#now that it is sorted in reverse put in bins
countY=1
for binY in orderedY
currentSize=binY[1].widths[1]
currentSizeLayer=1
append!(orderedYandSize[countY],[[]]);
for cube in binY
if(currentSize==cube.widths[1])
append!(orderedYandSize[countY][currentSizeLayer],[cube]);
else
append!(orderedYandSize[countY],[[]]);
currentSize=cube.widths[1]
currentSizeLayer=currentSizeLayer+1
append!(orderedYandSize[countY][currentSizeLayer],[cube]);
end
end
countY+=1
end
#for each orderY bin and size bin order according to distance value
orderedYandSizeandSDF=[]
countY=1
for binY in orderedYandSize
append!(orderedYandSizeandSDF,[[]]);
countSize=1
for binSize in binY
append!(orderedYandSizeandSDF[countY],[[]]);
orderedYandSize[countY][countSize]=sort(binSize,by = x -> getSDF(x))
countSize+=1
end
countY+=1