Commit 8012b5fe authored by Amira Abdel-Rahman's avatar Amira Abdel-Rahman
Browse files

2d slices 3d concurrent

parent 3e3867e6
Pipeline #12137 passed with stage
in 34 seconds
This diff is collapsed.
This diff is collapsed.
......@@ -13645,7 +13645,7 @@
},
{
"cell_type": "code",
"execution_count": 371,
"execution_count": 375,
"metadata": {},
"outputs": [
{
......@@ -13654,7 +13654,7 @@
"GLMakie.Screen(...)"
]
},
"execution_count": 371,
"execution_count": 375,
"metadata": {},
"output_type": "execute_result"
}
......@@ -13663,9 +13663,9 @@
"Micro_xPhys=Micro_xPhys1\n",
"AbstractPlotting.inline!(false)\n",
"\n",
"n=1\n",
"n=3\n",
"# for i=1:θ\n",
"i=3\n",
"i=1\n",
"temp=copy(Micro_xPhys[i])\n",
"# display(sum(Micro_xPhys[i]))\n",
"temp[Micro_xPhys[i].<0.7].=0.0\n",
This diff is collapsed.
This diff is collapsed.
This source diff could not be displayed because it is too large. You can view the blob instead.
This source diff could not be displayed because it is too large. You can view the blob instead.
This diff is collapsed.
{"nodes":[{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[6,7,8,9,10,11],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n0","position":{"x":0,"z":5,"y":0},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":0,"restrained_degrees_of_freedom":[false,true,false,true,true,true],"displacement":{"x":-5639.776228563598,"z":2048.3814342322685,"y":0.0}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[60,61,62,63,64,65],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n1","position":{"x":5,"z":5,"y":0},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":1,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":-5639.800959557511,"z":955.1576294722763,"y":-1908.3414608424196}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[30,31,32,33,34,35],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n2","position":{"x":0,"z":10,"y":5},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":2,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":-493.84236552822745,"z":-6057.114081138011,"y":8105.525252654581}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[84,85,86,87,88,89],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n3","position":{"x":5,"z":10,"y":5},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":3,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":-493.83597558188603,"z":-3097.285568939793,"y":-0.15311890169900177}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[18,19,20,21,22,23],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n4","position":{"x":0,"z":0,"y":5},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":4,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":-3345.253368062673,"z":7316.414879779308,"y":5268.703530264545}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[36,37,38,39,40,41],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n5","position":{"x":0,"z":0,"y":10},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":5,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":9615.02290685203,"z":-748.9258706432531,"y":5268.718346081642}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[48,49,50,51,52,53],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n6","position":{"x":0,"z":10,"y":10},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":6,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":12834.53875792117,"z":-734.706513313299,"y":8105.535374002722}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[66,67,68,69,70,71],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n7","position":{"x":5,"z":10,"y":0},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":7,"restrained_degrees_of_freedom":[false,true,false,true,true,true],"displacement":{"x":-8599.385479659079,"z":85.50033739033768,"y":0.0}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[42,43,44,45,46,47],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n8","position":{"x":0,"z":5,"y":10},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":1,"z":0,"y":-1},"idd":8,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":13471.296907985694,"z":-923.952294358437,"y":13238.678978684007}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[102,103,104,105,106,107],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n9","position":{"x":5,"z":10,"y":10},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":9,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":14406.659210085181,"z":-1859.3109483151886,"y":1055.7209052709586}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[54,55,56,57,58,59],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n10","position":{"x":5,"z":0,"y":0},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":10,"restrained_degrees_of_freedom":[false,true,false,true,true,true],"displacement":{"x":-7735.168375385161,"z":1238.0795337388813,"y":0.0}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[72,73,74,75,76,77],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n11","position":{"x":5,"z":0,"y":5},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":11,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":-5347.618030608182,"z":2340.353544855823,"y":-0.15777216904301652}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[120,121,122,123,124,125],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n12","position":{"x":10,"z":10,"y":0},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":12,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":-8599.398447489326,"z":3914.7223890498412,"y":-8105.772309129612}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[144,145,146,147,148,149],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n13","position":{"x":10,"z":0,"y":10},"force":{"x":0,"z":0,"y":-10},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":13,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":19522.524313395068,"z":2170.216673491202,"y":-24870.345125154476}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[78,79,80,81,82,83],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n14","position":{"x":5,"z":5,"y":5},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":14,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":381.38429984858584,"z":60.46606280795277,"y":381.29522790906515}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[114,115,116,117,118,119],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n15","position":{"x":10,"z":5,"y":0},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":15,"restrained_degrees_of_freedom":[true,true,true,true,true,true],"displacement":{"x":0.0,"z":0.0,"y":0.0}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[138,139,140,141,142,143],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n16","position":{"x":10,"z":10,"y":5},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":16,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":-499.8492528101398,"z":941.7093151141274,"y":-8105.898922745866}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[156,157,158,159,160,161],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n17","position":{"x":10,"z":10,"y":10},"force":{"x":0,"z":0,"y":-10},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":17,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":24310.473084096622,"z":-191.2647988403951,"y":-24804.66776980199}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[90,91,92,93,94,95],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n18","position":{"x":5,"z":0,"y":10},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":18,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":12037.384017452787,"z":-1859.2957180133815,"y":-3859.570417232936}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[96,97,98,99,100,101],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n19","position":{"x":5,"z":5,"y":10},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":19,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":15492.914148645705,"z":-1859.3212708394676,"y":868.9456797686016}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[24,25,26,27,28,29],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n20","position":{"x":0,"z":5,"y":5},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":20,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":1672.15212569106,"z":158.6335954867011,"y":8523.521275884286}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[132,133,134,135,136,137],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n21","position":{"x":10,"z":5,"y":5},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":21,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":2017.1835764689276,"z":-586.2799797741986,"y":-12606.764734374341}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[108,109,110,111,112,113],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n22","position":{"x":10,"z":0,"y":0},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":22,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":-6594.9883719792615,"z":-0.010713925673371705,"y":-13193.037141124321}},{"parent":"11","nomSize":1,"angle":{"x":0,"z":0,"y":0},"degrees_of_freedom":[126,127,128,129,130,131],"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"id":"n23","position":{"x":10,"z":0,"y":5},"force":{"x":0,"z":0,"y":0},"in":true,"fixedDisplacement":{"x":0,"z":0,"y":0},"idd":23,"restrained_degrees_of_freedom":[false,false,false,true,true,true],"displacement":{"x":-2168.446746342561,"z":-2489.372484384959,"y":-13193.05196017522}}],"hierarchical":false,"animation":{"exaggeration":2000,"speed":3,"showDisplacement":true},"viz":{"colorMap":0,"colorMaps":[],"maxStress":-10000000,"minStress":10000000},"edges":[{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":0,"id":"e8","stress":0,"target":1,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":0,"id":"e11","stress":0,"target":2,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":0,"id":"e13","stress":0,"target":3,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":4,"id":"e17","stress":0,"target":5,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":4,"id":"e24","stress":0,"target":0,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":2,"id":"e37","stress":0,"target":6,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":2,"id":"e38","stress":0,"target":3,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":2,"id":"e40","stress":0,"target":7,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":8,"id":"e50","stress":0,"target":9,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":8,"id":"e51","stress":0,"target":2,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":10,"id":"e55","stress":0,"target":11,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":1,"id":"e69","stress":0,"target":12,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":7,"id":"e73","stress":0,"target":3,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":7,"id":"e74","stress":0,"target":12,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":11,"id":"e78","stress":0,"target":13,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":11,"id":"e86","stress":0,"target":0,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":14,"id":"e92","stress":0,"target":15,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":14,"id":"e95","stress":0,"target":16,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":3,"id":"e104","stress":0,"target":17,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":3,"id":"e105","stress":0,"target":12,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":18,"id":"e108","stress":0,"target":19,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":18,"id":"e112","stress":0,"target":20,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":18,"id":"e113","stress":0,"target":21,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":19,"id":"e115","stress":0,"target":21,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":19,"id":"e116","stress":0,"target":9,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":22,"id":"e124","stress":0,"target":23,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":22,"id":"e125","stress":0,"target":15,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":22,"id":"e126","stress":0,"target":21,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":22,"id":"e127","stress":0,"target":1,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":12,"id":"e134","stress":0,"target":16,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":23,"id":"e139","stress":0,"target":14,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":21,"id":"e145","stress":0,"target":12,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":21,"id":"e146","stress":0,"target":3,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":21,"id":"e148","stress":0,"target":9,"targetNodalCoordinate":0},{"sourceNodalCoordinate":0,"material":{"density":0.028,"area":1.1328799999999999,"stiffness":4000},"source":13,"id":"e152","stress":0,"target":19,"targetNodalCoordinate":0}],"ndofs":162}
\ No newline at end of file
This diff is collapsed.
......@@ -22,6 +22,14 @@
# dispFlag = 1;
# plotFlag = 1;
# props, SH = evaluateCH(CH, dens, outOption, dispFlag);
# using MAT
using LinearAlgebra
# using CUDA, IterativeSolvers
using SparseArrays
using Krylov
using VectorizedRoutines
using GLMakie
function KE_from_E(E,nu)
D0 = E/(1+nu)/(1-2*nu)*
[ 1-nu nu nu 0 0 0 ;
......@@ -429,6 +437,200 @@ function evaluateCH(CH, dens, outOption, dispFlag)
return props, SH
end
###############################################
# input the homogenized elasticity tensor
function visual(CH)
# transform it to 3*3*3*3 tensor
tensor = generate(CH);
# find the E1 in 360 degree
# x = 0:pi/180:2*pi;
a,e = Matlab.meshgrid(0:0.02*π:2*π, -π/2:0.01*π:π/2);
E1 = zeros(size(a)[1],size(a)[2]);
for i = 1:size(a)[1]
for j = 1:size(a)[2]
# build the transformation matrix
trans_z = [cos(a[i,j]) -sin(a[i,j]) 0;
sin(a[i,j]) cos(a[i,j]) 0;
0 0 1];
trans_y = [cos(e[i,j]) 0 sin(e[i,j]);
0 1 0;
-sin(e[i,j]) 0 cos(e[i,j])];
# calculate the new tensor
N_tensor = transform(tensor,trans_y*trans_z);
# transform the tensor to 6*6
N_CH = ToMatrix(N_tensor);
# calculate the E1
E = modulus(N_CH);
E1[i,j] = E[1];
end
end
x,y,z = sph2cart(a,e,E1);
c = sqrt.(x.^2 .+ y.^2 +z.^2)
display(surface(x,y,z ,color=c))
return s
end
function modulus(CH)
S = inv(CH);
E = zeros(6,1);
E[1] = 1/S[1,1];
E[2] = 1/S[2,2];
E[3] = 1/S[3,3];
E[4] = 1/S[4,4];
E[5] = 1/S[5,5];
E[6] = 1/S[6,6];
return E
end
function generate(CH)
C = zeros(3,3,3,3);
for i = 1:6
for j = 1:6
a,b = change(i);
c,d = change(j);
C[a,b,c,d] = CH[i,j];
end
end
for i = 1:3
if (i == 3)
j = 1;
else
j = i+1;
end
for m = 1:3
if (m == 3)
n = 1;
else
n = m+1;
end
C[j,i,n,m] = C[i,j,m,n];
C[j,i,m,n] = C[i,j,m,n];
C[i,j,n,m] = C[i,j,m,n];
C[j,i,m,m] = C[i,j,m,m];
C[m,m,j,i] = C[m,m,i,j];
end
end
return C
end
# change the index 4 5 6 to 23 31 12
function change(w)
if (w < 4)
a = w;
b = w;
else
if (w == 4)
a = 2;
b = 3;
else
if (w == 5)
a = 3;
b = 1;
else
if (w==6)
a = 1;
b = 2;
end
end
end
end
return a,b
end
function ToMatrix(C)
CH2 = zeros(6,6);
for i = 1:6
for j = 1:6
a,b = change(i);
c,d = change(j);
CH2[i,j] = C[a,b,c,d];
end
end
return CH2
end
function transform(itr,tmx)
#
# FUNCTION
# otr = transform(itr,tmx)
#
# DESCRIPTION
# transform 3D-tensor (Euclidean or Cartesion tensor) of any order (>0) to another coordinate system
#
# PARAMETERS
# otr = output tensor, after transformation; has the same dimensions as the input tensor
# itr = input tensor, before transformation; should be a 3-element vector, a 3x3 matrix, or a 3x3x3x... multidimensional array, each dimension containing 3 elements
# tmx = transformation matrix, 3x3 matrix that contains the direction cosines between the old and the new coordinate system
#
ne = length(itr); # number of tensor elements
nd = ndims(itr); # number of tensor dimensions, i.e. order of tensor
if ne==3
nd = 1;
end # order of tensor is 1 in case of a 3x1 or 1x3 vector
otr = copy(itr); # create output tensor
otr[:] .= 0; # fill output tensor with zeros; this way a symbolic tensor remains symbolic
iie = zeros(nd); # initialise vector with indices of input tensor element
ioe = zeros(nd); # initialise vector with indices of output tensor element
cne = Int.(cumprod(3.0*ones(nd),dims=1)./3.0); # vector with cumulative number of elements for each dimension (divided by three)
# display(cne)
for oe = 1:ne
ioe =(floor.((oe .-1)./cne) .%3).+1; # calculate indices of current output tensor element
for ie = 1:ne # loop over all input elements
pmx = 1; # initialise product of transformation matrices
iie = (floor.((ie-1)./cne) .%3).+1; # calculate indices of current input tensor element
for id = 1:nd # loop over all dimensions
pmx = pmx * tmx[Int(ioe[id]), Int(iie[id])] ;
end
otr[oe] = otr[oe] + pmx * itr[ie]; # add product of transformation matrices and input tensor element to output tensor element
end
end
# Transform matrix about Z axis
# for x = 0:pi/20:pi/2
# trans = [cos(x) cos(x-pi/2) 0;
# cos(x+pi/2) cos(x) 0;
# 0 0 1];
# N_CH = transform(C,trans);
# end
return otr
end
function sph2cart(az,elev,r)
#SPH2CART Transform spherical to Cartesian coordinates.
# [X,Y,Z] = SPH2CART(TH,PHI,R) transforms corresponding elements of
# data stored in spherical coordinates (azimuth TH, elevation PHI,
# radius R) to Cartesian coordinates X,Y,Z. The arrays TH, PHI, and
# R must be the same size (or any of them can be scalar). TH and
# PHI must be in radians.
#
# TH is the counterclockwise angle in the xy plane measured from the
# positive x axis. PHI is the elevation angle from the xy plane.
#
# Class support for inputs TH,PHI,R:
# float: double, single
#
# See also CART2SPH, CART2POL, POL2CART.
# Copyright 1984-2005 The MathWorks, Inc.
z = r .* sin.(round.(elev,digits=16));
rcoselev = r .* cos.(round.(elev,digits=16));
x = rcoselev .* cos.(round.(az,digits=16));
y = rcoselev .* sin.(round.(az,digits=16));
return x,y,z
end
###############################################
## SUB FUNCTION: elementMatVec3D
function elementMatVec3D(a, b, c, DH)
......
......@@ -55,8 +55,11 @@ include("./problems.jl")
include("./concurrent2D.jl")
include("./concurrent3D.jl")
include("./fabrication.jl")
include("./topmgcg.jl")
println("Loaded Topology Optimization Library!")
#======================================================================================================================#
# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2021
#======================================================================================================================#
# using Clustering
# using StatsPlots
......@@ -399,8 +401,8 @@ function MultiConTop2DU(θ,Macro_struct, Micro_struct,prob, penal,saveItr=5,maxl
edofVec = reshape(2*nodenrs[1:end-1,1:end-1].+1,Macro_nele,1);
edofMat = repeat(edofVec,1,8).+repeat([0 1 2*Macro_nely.+[2 3 0 1] -2 -1],Macro_nele,1)
Micro_rmin=Int(Micro_struct[6])
Macro_rmin=Int(Macro_struct[6])
Micro_rmin=(Micro_struct[6])
Macro_rmin=(Macro_struct[6])
U,F,freedofs=prob(Macro_nelx,Macro_nely)
......@@ -408,8 +410,8 @@ function MultiConTop2DU(θ,Macro_struct, Micro_struct,prob, penal,saveItr=5,maxl
iK = reshape(kron(edofMat,ones(8,1))',64*Macro_nele,1);
jK = reshape(kron(edofMat,ones(1,8))',64*Macro_nele,1);
# PREPARE FILTER
Macro_H,Macro_Hs = filtering2d(Macro_nelx, Macro_nely, Macro_nele, Macro_rmin);
Micro_H,Micro_Hs = filtering2d(Micro_nelx, Micro_nely, Micro_nele, Micro_rmin);
Macro_H,Macro_Hs = make_filter(Macro_nelx, Macro_nely, Macro_rmin);
Micro_H,Micro_Hs = make_filter(Micro_nelx, Micro_nely, Micro_rmin);
## Connectivity matrix for continuity
L=connectivityMatrix(Micro_nelx,Micro_nely);
......@@ -494,7 +496,8 @@ function MultiConTop2DU(θ,Macro_struct, Micro_struct,prob, penal,saveItr=5,maxl
end
loopbeta = 0; loop = 0; Macro_change = 1;
while loop < maxloop && ( Macro_change > 0.01 || Bool.(sum(Micro_change .>= 0.01)<θ))
while loop < maxloop && ( Macro_change > 0.01 || sum(Micro_change .>= 0.01)>=1)
# while loop < maxloop || Macro_change > 0.01 || Bool.(sum(Micro_change .>= 0.01))
loop = loop+1; loopbeta = loopbeta+1;
......@@ -578,7 +581,7 @@ function MultiConTop2DU(θ,Macro_struct, Micro_struct,prob, penal,saveItr=5,maxl
ss=6
Micro_x[i]=addFabConstraint2D(Micro_x[i],ss);
Micro_xPhys[i]=addFabConstraint2D(Micro_xPhys[i],ss);
Micro_x[i], Micro_xPhys[i], Micro_change1 = OC2D(Micro_x[i], Micro_dc[i], Micro_dv[i], Micro_H, Micro_Hs, Micro_Vol[i], Micro_nele, 0.2, beta);
Micro_x[i], Micro_xPhys[i], Micro_change1 = OC2D(Micro_x[i], Micro_dc[i], Micro_dv[i], Micro_H, Micro_Hs, Micro_Vol[i], Micro_nele, 0.2, beta,true,ss);
# Micro_x[i], Micro_xPhys[i], Micro_change1,Micro_low[i],Micro_upp[i],Micro_xold1[i],Micro_xold2[i]= OC2DMMA(Micro_x[i],Micro_xold1[i],Micro_xold2[i], Micro_dc[i], Micro_dv[i], Micro_H, Micro_Hs, Micro_Vol[i], Micro_nele, 0.2, beta,loop,Micro_low[i],Micro_upp[i]);
else
Micro_x[i], Micro_xPhys[i], Micro_change1,Micro_low[i],Micro_upp[i],Micro_xold1[i],Micro_xold2[i]= OC2DMMACON(Micro_x[i],Micro_xold1[i],Micro_xold2[i], Micro_dc[i], Micro_dv[i], Micro_H, Micro_Hs, Micro_Vol[i], Micro_nele, 0.2, beta,loop,Micro_low[i],Micro_upp[i],L);
......@@ -781,7 +784,7 @@ function CompliantMultiConTop2DU(θ,Macro_struct, Micro_struct,prob, penal,saveI
# loopbeta = 0; loop = 0; Macro_change = 1;
for penal=1.0:4
loopbeta = 0; loop = 0; Macro_change = 1;
while loop < maxloop && ( Macro_change > 0.01 || Bool.(sum(Micro_change .>= 0.01)<θ))
while loop < maxloop && ( Macro_change > 0.01 || sum(Micro_change .>= 0.01)>=1)
loop = loop+1; loopbeta = loopbeta+1;
# FE-ANALYSIS AT TWO SCALES
......@@ -1074,11 +1077,131 @@ function Uclustering2D(θ,Macro_nely,Macro_nelx,Macro_nele,Macro_Vol,Micro_Vol,i
K = sparse(vec(iK),vec(jK),vec(sK)); K = (K+K')/2;
U[freedofs,:] = K[freedofs,freedofs]\Array(F[freedofs,:]);
########################################################################
#visualize deformation/strain
# display(Plots.heatmap(1.0.-Macro_xPhys[end:-1:1,:],showaxis = false,xaxis=nothing,yaxis=nothing,legend=nothing,clims=(0, 1),fc=:grays,aspect_ratio=:equal))
diss=reshape(U[edofMat],Macro_nely, Macro_nelx,8);order=[1,2,3,4,5,6,7,8];
# gr(size=(50*4,50*3))
# for j=1:8
# display(Plots.heatmap(diss[:,:,j],showaxis = false,xaxis=nothing,yaxis=nothing,legend = false,aspect_ratio=:equal))
# end
# gr(size=(400,300))
unit=10
ny=size(diss)[1];nx=size(diss)[2];
ϵx =zeros(ny ,nx);ϵy =zeros(ny ,nx);
αx =zeros(ny ,nx);αy =zeros(ny ,nx);
α1 =zeros(ny ,nx);α2 =zeros(ny ,nx);α3 =zeros(ny ,nx);α4 =zeros(ny ,nx);
ϵxy =zeros(ny ,nx);τxy=zeros(ny ,nx);
u=zeros(ny ,nx);v=zeros(ny ,nx);
Θ=zeros(ny ,nx);ΘD=zeros(ny ,nx);
von_mises=zeros(ny ,nx);
l=100;nu=0.3;ν=0.3;E=100;
a=1;b=1;
x=0;y=0;#at center
#shape functions
N1=(a-x)*(b-y)/(4*a*b)
N2=(a+x)*(b-y)/(4*a*b)
N3=(a+x)*(b+y)/(4*a*b)
N4=(a-x)*(b+y)/(4*a*b)
B = (1/2/l)*[-1 0 1 0 1 0 -1 0; 0 -1 0 -1 0 1 0 1; -1 -1 -1 1 1 1 1 -1];
C = E/(1-ν^2)*[1 ν 0; ν 1 0; 0 0 (1-ν)/2];
C = E/((1+ν)*(1-2*ν))*[1-ν ν 0; ν 1-ν 0; 0 0 (1-2*ν)/2];
for i=1:ny
for j=1:nx
#based on: https://community.wvu.edu/~bpbettig/MAE456/Exam_Final_practice_answers.pdf
Ue = diss[i,j,:];
N=[N1 0 N2 0 N3 0 N4 0; 0 N1 0 N2 0 N3 0 N4];
#displacement at center
u[i,j],v[i,j]=N*Ue
#in plane strain in center
ϵ=1/(4*a*b)*[-(b-y) 0 (b-y) 0 (b+y) 0 -(b+y) 0;
0 -(a-x) 0 -(a+x) 0 (a+x) 0 (a-x);
-(a-x) -(b-y) -(a+x) (b-y) (a+x) (b+y) (a-x) -(b+y)]*Ue
ϵx[i,j],ϵy[i,j],ϵxy[i,j]=ϵ
#in plane stress in center
α=C*ϵ
αx[i,j],αy[i,j],τxy[i,j]=α
#principal stress center
α1[i,j]=(αx[i,j]+αy[i,j])/2+sqrt(((αx[i,j]-αy[i,j])/2)^2+τxy[i,j]^2)
α2[i,j]=(αx[i,j]+αy[i,j])/2-sqrt(((αx[i,j]-αy[i,j])/2)^2+τxy[i,j]^2)
α3[i,j]=ν*(αx[i,j]+αy[i,j])
α4[i,j]=ν*(αx[i,j]+αy[i,j])
# Θ =atan((2*τxy),(αx-αy))/2
# Θ =atan((2*ϵxy),(ϵx-ϵy))/2
#von_mises
von_mises[i,j]=1/(sqrt(2))*sqrt((α1[i,j]-α2[i,j])^2+(α2[i,j]-α3[i,j])^2+(α3[i,j]-α1[i,j])^2)
Θ[i,j] =atan((2*ϵxy[i,j]/2),(ϵx[i,j]-ϵy[i,j]))/2
Θ[i,j] =atan((2*τxy[i,j]),(αx[i,j]-αy[i,j]))/2
if abs(α1[i,j])<abs(α2[i,j])
Θ[i,j] -=π/2
end
end
end
ΘD=Θ.*180/π;
Plots.heatmap(von_mises, aspect_ratio=:equal)#,clim=(0,2e-5))
xs = 1:Int(nx)
ys = 1:Int(ny)
maxi=maximum(abs.(α1))/2
# ress(y,x)=(cos(Θ[Int(x),Int(y)])*α1[Int(x),Int(y)]/maxi, sin(Θ[Int(x),Int(y)])*α1[Int(x),Int(y)]/maxi)
# ress1(y,x)=(cos(Θ[Int(x),Int(y)]-π/2)*α2[Int(x),Int(y)]/maxi, sin(Θ[Int(x),Int(y)]-π/2)*α2[Int(x),Int(y)]/maxi)
ress(y,x)=(cos(Θ[Int(x),Int(y)])/3, sin(Θ[Int(x),Int(y)])/3)
ress1(y,x)=(cos(Θ[Int(x),Int(y)]+π/2)/3, sin(Θ[Int(x),Int(y)]+π/2)/3)
ress2(y,x)=(cos(Θ[Int(x),Int(y)]-π/2)/3, sin(Θ[Int(x),Int(y)]-π/2)/3)
ress3(y,x)=(cos(Θ[Int(x),Int(y)]+π)/3, sin(Θ[Int(x),Int(y)]+π)/3)
# col(y,x)=norm(ress(y,x))
xxs = [x for x in xs for y in ys]
yys = [y for x in xs for y in ys]
# Plots.quiver!(xxs, yys, quiver=ress)
# Plots.quiver!(xxs, yys, quiver=ress1)
# Plots.quiver!(xxs, yys, quiver=ress2)
# display(Plots.quiver!(xxs, yys, quiver=ress3))
display(Plots.quiver!(xxs, yys, quiver=ress))
dissVec=zeros(ny+1 ,nx+1 ,2)
for i=1:ny
for j=1:nx
dissVec[i,j,1]=diss[i,j,1]
dissVec[i,j,2]=diss[i,j,2]
end
end
for i=1:ny
dissVec[i,end,1]=diss[i,end,order[3]]
dissVec[i,end,2]=diss[i,end,order[4]]
end
for j=1:nx
dissVec[end,j,1]=diss[end,j,order[7]]
dissVec[end,j,2]=diss[end,j,order[8]]
end
dissVec[end,end,1]=diss[end,end,order[5]]
dissVec[end,end,2]=diss[end,end,order[6]]
UX=dissVec[:,:,1];UY=dissVec[:,:,2]
plotStrain(nx,ny,unit,diss,order,UX,UY)
########################################################################
# X=U[edofMat]'
X=hcat(Macro_xPhys[:],U[edofMat])'
# X=hcat(Macro_xPhys[:],U[edofMat])'
X=hcat(Macro_xPhys[:],Θ[:])'
# X=hcat(Θ[:])'
# cluster X into θ clusters using K-means
# R = kmeans(X, θ; maxiter=200, display=:iter)
R = Clustering.kmeans(X, θ; maxiter=200)
R = Clustering.kmeans(X, θ; maxiter=1000)
@assert Clustering.nclusters(R) == θ # verify the number of clusters
......@@ -1088,7 +1211,7 @@ function Uclustering2D(θ,Macro_nely,Macro_nelx,Macro_nele,Macro_Vol,Micro_Vol,i
D=Distances.pairwise(Euclidean(), M, dims=2)
hc = hclust(D, linkage=:single)
display(Plots.plot(hc))
Uclustered=reshape(a,Macro_nely,Macro_nelx);
# display(Plots.heatmap(Uclustered,showaxis = false,xaxis=nothing,yaxis=nothing,clims=(0.0, θ),fc=jet_me,aspect_ratio=:equal))
......@@ -1124,7 +1247,7 @@ function Uclustering2D(θ,Macro_nely,Macro_nelx,Macro_nele,Macro_Vol,Micro_Vol,i
# display(Macro_masks_dens)
display(Plots.heatmap(Macro_xPhys_diagU,showaxis = false,xaxis=nothing,yaxis=nothing,clims=(1, θ),fc=cscheme,aspect_ratio=:equal))
display(Plots.plot(hc))
return Macro_masksU,Macro_masks_densU,Macro_xPhys_diagU
end
......@@ -1193,6 +1316,123 @@ function CompliantUclustering2D(θ,Macro_nely,Macro_nelx,Macro_nele,Macro_Vol,Mi
sK = reshape(Ke[:]*(Emin .+Macro_xPhys[:]'.^penal1*(1 .-Emin)),64*Macro_nele,1);
K = sparse(vec(iK),vec(jK),vec(sK)); K = (K+K')/2;
U[freedofs,:] = K[freedofs,freedofs]\Array(F[freedofs,:]);
########################################################################
#visualize deformation/strain
# display(Plots.heatmap(1.0.-Macro_xPhys[end:-1:1,:],showaxis = false,xaxis=nothing,yaxis=nothing,legend=nothing,clims=(0, 1),fc=:grays,aspect_ratio=:equal))