// Create Sphero Cylindrical Fibers in GMSH using OpenCASCADE! SetFactory("OpenCASCADE"); // Set mesh charateristics for mesh not on fiber surface //h = 0.5; //Mesh.CharacteristicLengthMin = h; //Mesh.CharacteristicLengthMax = h; // Ask OpenCASCADE to compute more accurate bounding boxes of entities using the // STL mesh: RefineMesh; Mesh.SecondOrderLinear=1; Geometry.OCCBoundsUseStl = 1; Geometry.ToleranceBoolean = 1e-20; // des te lager des te beter, 1 gaat meteenfout e-30 lijkt het beste Geometry.MatchGeomAndMesh = 1; // Dit heeft positieve invloed //Geometry.MatchMeshTolerance = 1e-20; //Geometry.Tolerance = 1e-4; Geometry.OCCAutoFix = 0; Geometry.OCCFixSmallEdges = 1; Geometry.OCCFixSmallFaces = 1; eps = 1e-3; allSurfsPrev = 0; ///////////// Functions /////////////// // Make Sphero-Cylindric Fiber Volume Function MakeSpheroCylindricFbrVolume Cylinder(1) = {x-(0.5*n1x*hr - 2^0.5*r*n1x),y-(0.5*n1y*hr - 2^0.5*r*n1y),z-(0.5*n1z*hr - 2^0.5*r*n1z), n1x*hr-2*2^0.5*r*n1x, n1y*hr-2*2^0.5*r*n1y, n1z*hr-2*2^0.5*r*n1z, r}; Sphere(2) = {x-(0.5*n1x*hr - 2^0.5*r*n1x),y-(0.5*n1y*hr - 2^0.5*r*n1y),z-(0.5*n1z*hr - 2^0.5*r*n1z), r}; Sphere(3) = {x+(0.5*n1x*hr - 2^0.5*r*n1x),y+(0.5*n1y*hr - 2^0.5*r*n1y),z+(0.5*n1z*hr - 2^0.5*r*n1z), r}; //100+3*(t-1)+ // add all volumes to create the sphero cylinder // All volumes involved are being deleted. BooleanUnion( boxIdx + t ) = { Volume{1}; Delete; }{ Volume{2,3}; Delete; }; allSurfs() = Surface In BoundingBox{xmin-2, ymin-2, zmin-2, xmax+4, ymax+4, zmax+4}; Printf("XXXXXtotalnr of surfaces '%g',", #allSurfs()); Return // Get Boundary Surfaces Function GetBoundarySurfaces Printf("XXXXXXpunt4 '%g','%g','%g',", xp(4), yp(4), zp(4)); Printf("XXXXXXpunt7 '%g','%g','%g',", xp(7), yp(7), zp(7)); For i In {0:#allSurfs()-1} // draw a box around the surface, if the box coordinates correspond to the corner points // of the domain then that boundary surface is found. bbox() = BoundingBox Surface { allSurfs(i) }; Printf("XXXXXXbb1 '%g',", bbox(0)); Printf("XXXXXXbb2 '%g',", bbox(1)); Printf("XXXXXXbb3 '%g',", bbox(2)); Printf("XXXXXXbb4 '%g',", bbox(3)); Printf("XXXXXXbb5 '%g',", bbox(4)); Printf("XXXXXXbb6 '%g',", bbox(5)); If (bbox(0) <= xp(1)+eps && bbox(0) >= xp(1)-eps && bbox(1) <= yp(1)+eps && bbox(1) >= yp(1)-eps && bbox(2) <= zp(1)+eps && bbox(2) >= zp(1)-eps && bbox(3) <= xp(3)+eps && bbox(3) >= xp(3)-eps && bbox(4) <= yp(3)+eps && bbox(4) >= yp(3)-eps && bbox(5) <= zp(3)+eps && bbox(5) >= zp(3)-eps ) S1(0) = allSurfs(i); Printf("XXXXXXS1 '%g',", S1(0)); EndIf If (bbox(0) <= xp(2)+eps && bbox(0) >= xp(2)-eps && bbox(1) <= yp(2)+eps && bbox(1) >= yp(2)-eps && bbox(2) <= zp(2)+eps && bbox(2) >= zp(2)-eps && bbox(3) <= xp(7)+eps && bbox(3) >= xp(7)-eps && bbox(4) <= yp(7)+eps && bbox(4) >= yp(7)-eps && bbox(5) <= zp(7)+eps && bbox(5) >= zp(7)-eps ) S2(0) = allSurfs(i); Printf("XXXXXXS2 '%g',", S2(0)); EndIf If (bbox(0) <= xp(5)+eps && bbox(0) >= xp(5)-eps && bbox(1) <= yp(5)+eps && bbox(1) >= yp(5)-eps && bbox(2) <= zp(5)+eps && bbox(2) >= zp(5)-eps && bbox(3) <= xp(7)+eps && bbox(3) >= xp(7)-eps && bbox(4) <= yp(7)+eps && bbox(4) >= yp(7)-eps && bbox(5) <= zp(7)+eps && bbox(5) >= zp(7)-eps ) S3(0) = allSurfs(i); Printf("XXXXXXS3 '%g',", S3(0)); EndIf If (bbox(0) <= xp(1)+eps && bbox(0) >= xp(1)-eps && bbox(1) <= yp(1)+eps && bbox(1) >= yp(1)-eps && bbox(2) <= zp(1)+eps && bbox(2) >= zp(1)-eps && bbox(3) <= xp(8)+eps && bbox(3) >= xp(8)-eps && bbox(4) <= yp(8)+eps && bbox(4) >= yp(8)-eps && bbox(5) <= zp(8)+eps && bbox(5) >= zp(8)-eps ) S4(0) = allSurfs(i); Printf("XXXXXXS4 '%g',", S4(0)); EndIf If (bbox(0) <= xp(4)+eps && bbox(0) >= xp(4)-eps && bbox(1) <= yp(4)+eps && bbox(1) >= yp(4)-eps && bbox(2) <= zp(4)+eps && bbox(2) >= zp(4)-eps && bbox(3) <= xp(7)+eps && bbox(3) >= xp(7)-eps && bbox(4) <= yp(7)+eps && bbox(4) >= yp(7)-eps && bbox(5) <= zp(7)+eps && bbox(5) >= zp(7)-eps ) S5(0) = allSurfs(i); Printf("XXXXXXS5 '%g',", S5(0)); EndIf If (bbox(0) <= xp(1)+eps && bbox(0) >= xp(1)-eps && bbox(1) <= yp(1)+eps && bbox(1) >= yp(1)-eps && bbox(2) <= zp(1)+eps && bbox(2) >= zp(1)-eps && bbox(3) <= xp(6)+eps && bbox(3) >= xp(6)-eps && bbox(4) <= yp(6)+eps && bbox(4) >= yp(6)-eps && bbox(5) <= zp(6)+eps && bbox(5) >= zp(6)-eps ) S6(0) = allSurfs(i); Printf("XXXXXXS6 '%g',", S6(0)); // make a Physical Point on this surface in order to set a essential boundary condition POS6() = PointsOf{ Surface{ allSurfs(i) }; }; //Physical Point(1) = POS6(0); EndIf EndFor Printf("XXXXXXS1 '%g',", S1(0)); Printf("XXXXXXS2 '%g',", S2(0)); Printf("XXXXXXS3 '%g',", S3(0)); Printf("XXXXXXS4 '%g',", S4(0)); Printf("XXXXXXS5 '%g',", S5(0)); Printf("XXXXXXS6 '%g',", S6(0)); Return //////////// End - Functions ////////// //------------------------------------------------------ //////////// Start of code ////////// // make the volume of the external domain, EVERYTHING IN OPENCASCADE SYNTAX! // In shear this is not a rectangular box anymore, we need to deal with this! boxIdx = 6; Box(1) = {xmin,ymin,zmin,xmax-xmin,ymax-ymin,zmax-zmin}; Box(2) = {xmin,ymin,zmin, xmin+Cos(roty)*lx, ymin+ly, zmin+10*zmax}; // rotate around y axis //Box(3) = {xmin,ymin,zmin, Cos(rotz)*lx, 4*ymax, lz}; // rotate around z axis //Box(4) = {xmin,ymin,zmin, lx, 4*ymax, Cos(rotx)*lz}; // rotate around x axis // Rotate additional boxes: Rotate {{0,1,0}, {xmin, ymin, zmin}, roty} { Volume{2}; } //Rotate {{0,0,1}, {xmin, ymin, zmin}, -rotz} { Volume{2}; } //Rotate {{1,0,0}, {xmin, ymin, zmin}, -roty} { Volume{4}; } // Intersection between boxes is the deformed box BooleanIntersection(boxIdx) = { Volume{2}; Delete; }{ Volume{1}; Delete; }; //BooleanIntersection(boxIdx) = { Volume{boxIdx-1}; Delete; }{ Volume{3,4}; Delete; }; // loop to make all spherocylidrical fibers For t In {1:nfibers} x = xr[t]; y = yr[t]; z = zr[t]; r = rfiber[t]*2^0.5; hr = hfiber[t]; n1x = n1rx[t]; n1y = n1ry[t]; n1z = n1rz[t]; Call MakeSpheroCylindricFbrVolume; EndFor a() = BooleanFragments{ Volume{boxIdx}; Delete; }{ Volume{boxIdx+1:boxIdx+nfibers}; Delete; }; Printf("XXXXXnew volumesl '%g',", #a()); For i In {0:#a()-1} Printf("XXXXXnew volumes '%g',", a(i)); // If (i!=0) // Delete { Volume{a(i)}; } // EndIf EndFor BooleanDifference( 100 ) = { Volume{a(0)}; Delete; }{ Volume{a(1):a(#a()-1)}; Delete; }; // Subtract all fibers volumes from the box volume //BooleanDifference( boxIdx+nfibers+1) = { Volume{boxIdx}; Delete; }{ Volume{boxIdx+1:boxIdx+nfibers}; Delete; }; //BooleanDifference( boxIdx+1) = { Volume{boxIdx}; Delete; }{ Volume{101:100+3*nfibers}; Delete; }; //Printf("XXXXXNU GAAN WE VOLUMES AFTREKKEN '%g',", #allSurfs()); //a() = BooleanDifference{ Volume{boxIdx}; Delete; }{ Volume{101:100+3*nfibers}; Delete; }; //Printf("XXXXXvolume a leng '%g',", #a()); //For t In {0:#a()-1} // Printf("XXXXXvolume a num '%g',", a(t)); //EndFor /// Next is finding the boundary surfaces to implement periodicity // get all surfaces in the domain allSurfs() = Surface In BoundingBox{xmin-2, ymin-2, zmin-2, xmax+4, ymax+4, zmax+4}; Printf("XXXXXtotalnr of surfaces_F '%g',", #allSurfs()); eps = 1; Call GetBoundarySurfaces; eps = 1e-3; // Implement periodicity //Periodic Surface {S2(0)} = {S4(0)} Translate {xp(2)-xp(1),yp(2)-yp(1),zp(2)-zp(1)}; //Periodic Surface {S3(0)} = {S1(0)} Translate {xp(5)-xp(1),yp(5)-yp(1),zp(5)-zp(1)}; //Periodic Surface {S5(0)} = {S6(0)} Translate {xp(4)-xp(1),yp(4)-yp(1),zp(4)-zp(1)}; // Create physical surfaces of the external domain //Physical Surface(1) = {S1(0)}; //Physical Surface(2) = {S2(0)}; //Physical Surface(3) = {S3(0)}; //Physical Surface(4) = {S4(0)}; //Physical Surface(5) = {S5(0)}; //Physical Surface(6) = {S6(0)}; // Get all fiber surfaces ( = all surfaces - external boundary surfaces) allFbrSurf() = allSurfs(); allFbrSurf() -=S1(0); allFbrSurf() -=S2(0); allFbrSurf() -=S3(0); allFbrSurf() -=S4(0); allFbrSurf() -=S5(0); allFbrSurf() -=S6(0); // At this point we have a set of randomly numbered surfaces and we need to connect // the correct surfaces to make a nice fiber. These are surfaces of the cylindre and // surfaces of the sphere. These surfaces can be connected because they share points! // SO find these matching points! numMatchingSurfaces = 0; matchesIdx = 0; For t In {0:#allFbrSurf()-1} POS1() = PointsOf{ Surface{ allFbrSurf(t) }; }; If (t < #allFbrSurf()-1) For k In {t+1:#allFbrSurf()-1} POS2() = PointsOf{ Surface{ allFbrSurf(k) }; }; For i In {0:#POS1()-1} For j In {0:#POS2()-1} If (POS1(i) == POS2(j)) // MATCHING POINT --> SURFACES ARE CONNECTED // CHECK IF CONNECTION IS ALREADY IN THE ARRAY If ( matchesIdx == 0 ) matches[matchesIdx] = allFbrSurf(t) ; matchesIdx = matchesIdx+1; matches[matchesIdx] = allFbrSurf(k) ; matchesIdx = matchesIdx+1; EndIf If( matches[matchesIdx-2] != allFbrSurf(t) || matches[matchesIdx-1] != allFbrSurf(k) ) matches[matchesIdx] = allFbrSurf(t) ; matchesIdx = matchesIdx+1; matches[matchesIdx] = allFbrSurf(k) ; matchesIdx = matchesIdx+1; EndIf EndIf EndFor EndFor EndFor EndIf EndFor // order and structurize fbrs = 0; surfStartArrIdx = 0; For k In {0:(#matches()/2)-1} If (k == 0) surfs(0) = matches(0); surfs(1) = matches(1); surfsIdx = 2; surfStartIdx = 0; continue = 1; EndIf If (k > 0 ) // CHECK IF matches(2*k) al in surfs zit continue = 1; For i In {0:#surfs()-1} If (matches(2*k) == surfs(i) ) continue = 0; EndIf EndFor If (continue == 1) surfs(surfsIdx) = matches(2*k); surfStartIdx = surfsIdx; surfsIdx = surfsIdx + 1; surfs(surfsIdx) = matches(2*k + 1 ); surfsIdx = surfsIdx + 1; EndIf EndIf If (continue == 1) surfStartArr[surfStartArrIdx] = surfStartIdx; surfStartArrIdx = surfStartArrIdx + 1; For i In {0:(#matches()/2)-1} If ( surfs(surfStartIdx) == matches(2*i) ) // already in surflist? insurflist = 0; For j In {surfStartIdx:#surfs()-1} // Aanpassen naar begin surface If ( surfs(j) == matches(2*i+1) ) // already in insurflist = 1; EndIf EndFor If (insurflist == 0) surfs(surfsIdx) = matches(2*i+1); surfsIdx = surfsIdx +1; EndIf EndIf EndFor EndIf // Loop over all added surfaces ( these can again match to a surface) // This should be a while loop/ or update the surfs() array For i In {surfStartIdx+1:#surfs()-1} For j In {0:(#matches()/2)-1} If ( surfs(i) == matches(2*j) ) // already in surflist? insurflist = 0; For ll In {surfStartIdx:#surfs()-1} // Aanpassen naar begin surface If ( surfs(ll) == matches(2*j+1) ) // already in insurflist = 1; EndIf EndFor If (insurflist == 0) surfs(surfsIdx) = matches(2*j+1); surfsIdx = surfsIdx +1; EndIf EndIf EndFor EndFor EndFor // Add the remaining fibers, which consist out of only one surface singleSurfFbrs() = allFbrSurf(); For i In {0:#surfs()-1} singleSurfFbrs() -= surfs(i); EndFor For i In {0:#singleSurfFbrs()-1} surfs(surfsIdx) = singleSurfFbrs(i); surfStartArr[surfStartArrIdx] = surfsIdx; surfsIdx = surfsIdx +1; surfStartArrIdx = surfStartArrIdx + 1; EndFor // length of the surfStartArr should be identical to the nfibers, check this! Printf("LENGTH OF THE SURFSTARTARR '%g',", #surfStartArr()); Printf("nfibers '%g',", nfibers); Printf("CHECK IF EQUAL!! '%g',", nfibers == #surfStartArr()); // Make the Physical surfaces and ensure they correspond to TFEM Fiber numbers For k In {0:#surfStartArr()-1} // Get point of surface, this point should be on the cylinder/spheros of the fibers POS() = PointsOf{ Surface{ surfs(surfStartArr(k)) }; }; /// LOOOP OVER AL DIE PUNTEN!!!! c() = Point{POS(0)}; Printf("XXXXXCOOR C: '%g','%g','%g',", c(0), c(1), c(2)); For t In {1:nfibers} x = xr[t]; y = yr[t]; z = zr[t]; r = rfiber[t]*2^0.5; hr = hfiber[t]; n1x = n1rx[t]; n1y = n1ry[t]; n1z = n1rz[t]; n2x = n2rx[t]; n2y = n2ry[t]; n2z = n2rz[t]; n3x = n3rx[t]; n3y = n3ry[t]; n3z = n3rz[t]; surfnr = surfacenr[t]; // Determine if point c is in the fiber(t) // fiber(t) consist out of cylinder with points p1c & p2c on the ends of the centerline // of the cylinder. p1c = { (x - (0.5*hr*n1x - 2^0.5*r*n1x)), (y - (0.5*hr*n1y - 2^0.5*r*n1y)), (z - (0.5*hr*n1z - 2^0.5*r*n1z))}; p2c = { (x + (0.5*hr*n1x - 2^0.5*r*n1x)), (y + (0.5*hr*n1y - 2^0.5*r*n1y)), (z + (0.5*hr*n1z - 2^0.5*r*n1z))}; p21c = { p2c(0)-p1c(0), p2c(1)-p1c(1), p2c(2)-p1c(2) }; // calculate the arm from centerline cylinder to the point c cfc_num(0) = ( c(1)-p1c(1) )*p21c(2) - ( c(2)-p1c(2) )*p21c(1); cfc_num(1) = ( c(2)-p1c(2) )*p21c(0) - ( c(0)-p1c(0) )*p21c(2); cfc_num(2) = ( c(0)-p1c(0) )*p21c(1) - ( c(1)-p1c(1) )*p21c(0); cfc_num_abs = (cfc_num(0)^2 + cfc_num(1)^2 + cfc_num(2)^2); cfc_den_abs = (p21c(0)^2 + p21c(1)^2 + p21c(2)^2); // point in top sphere, point in bottom sphere or point in cilinder If ( ( (c(0)-(x - (0.5*hr*n1x - 2^0.5*r*n1x)))^2 + (c(1)-(y - (0.5*hr*n1y - 2^0.5*r*n1y)) )^2 + (c(2)-(z - (0.5*hr*n1z - 2^0.5*r*n1z)))^2 < (r)^2+eps ) || ( (c(0)-(x + (0.5*hr*n1x - 2^0.5*r*n1x)))^2 + (c(1)-(y + (0.5*hr*n1y - 2^0.5*r*n1y)) )^2 + (c(2)-(z + (0.5*hr*n1z - 2^0.5*r*n1z)))^2 < (r)^2+eps ) || ( (c(0)-p1c(0))*p21c(0) + (c(1)-p1c(1))*p21c(1) + (c(2)-p1c(2))*p21c(2) >= 0 && (c(0)-p2c(0))*p21c(0) + (c(1)-p2c(1))*p21c(1) + (c(2)-p2c(2))*p21c(2) <= 0 && cfc_num_abs/cfc_den_abs <= (r)^2+eps ) ) Printf("link_TFEM FOUND: fiber '%g',", t ); Printf("Surface NR '%g',", surfnr ); If (k < #surfStartArr()-1) // Physical Surface(surfnr) = { surfs({surfStartArr(k):surfStartArr(k+1)-1}) }; EndIf If (k == #surfStartArr()-1) // Physical Surface(surfnr) = { surfs({ surfStartArr(k):#surfs()-1 }) }; EndIf EndIf EndFor EndFor // WHAT IF NO MATCH IS FOUND / two or more matches are found ( fiber fits in 2 surfaces??? ) // OUTPUT PRINTING FOR DEBUG For i In {0:(#matches()/2)-1} Printf("match '%g','%g',", matches(2*i), matches(2*i+1) ); EndFor Printf("surfslength '%g',", #surfs() ); For i In {0:#surfs()-1} Printf("ss1 '%g',", surfs(i)); EndFor Printf("ssALENGTH '%g',", #surfStartArr()); For i In {0:#surfStartArr()-1} Printf("ssA '%g',", surfStartArr(i)); EndFor // SET ELEMENT SIZE ON THE FIBER SURFACES! MeshSize{ PointsOf{ Surface { allFbrSurf() }; } } = 0.1; //Physical Volume(1) = {boxIdx+nfibers+1};