1
$\begingroup$
  Needs["TetGenLink`"];
    Needs["ComputationalGeometry`"];

    long = 100;
    broad = 100;
    height = 100;

    pts = {{0, 0, 0}, {long, 0, 0}, {long, broad, 0}, {0, broad, 0},
       {0, 0, height}, {long, 0, height}, {long, broad, height}, {0, 
        broad, height}};
    facets = {{{1, 2, 3, 4}}, {{5, 6, 7, 8}}, {{1, 5, 6, 2}}, {{2, 6, 7, 
         3}}, {{3, 7, 8, 4}}, {{4, 8, 5, 1}}};

    inst = TetGenCreate[];
    TetGenSetPoints[inst, pts];
    TetGenSetFacets[inst, facets];

    inst1 = TetGenTetrahedralize[inst, "pq1.414a4"];
    ElemPtsCoord = TetGenGetPoints[inst1];
    ElemFacesMadeofPtsIndex = TetGenGetFaces[inst1];
    ElemElemsMadeofPtsIndex = TetGenGetElements[inst1];

    Graphics3D[GraphicsComplex[ElemPtsCoord, Polygon[ElemFacesMadeofPtsIndex]]]

    ElemPtsCoordandIndex = ElemPtsCoord;
    Do[PrependTo[ElemPtsCoordandIndex[[i]], i], {i, 1, 
       Length[ElemPtsCoordandIndex]}
      ];

    ElemElemsMadeofPtsIndexandIdex = ElemElemsMadeofPtsIndex;
    Do[PrependTo[ElemElemsMadeofPtsIndexandIdex[[i]], i], {i, 1, 
       Length[ElemElemsMadeofPtsIndexandIdex]}
      ];

   RetainElementsIndex=Range[1,Length[ElemElemsMadeofPtsIndexandIdex]];

   ProjectedElemsIndex[ElemPtsCoord_, ElemElemsMadeofPtsIndexandIdex_, 
  Acen_, Bcen_, Ccen_, R1_, RetainElementsIndex_] := Module[
  {RetainElemsMadeofPtsIndexandIndex, ClassifiedElements, 
   ProjectedElementsIndex},

  RetainElemsMadeofPtsIndexandIndex = 
   ElemElemsMadeofPtsIndexandIdex[[RetainElementsIndex]];

  ClassifiedElements = GatherBy[
    RetainElemsMadeofPtsIndexandIndex,
    Volume[
       RegionIntersection[
        Ball[{Acen, Bcen, Ccen}, R1],
        Simplex[
         ElemPtsCoord[[  Delete[#1, 1]  ]]
         ]
        ]
       ] != 0 &];

  ProjectedElementsIndex = If[
    Volume[
      RegionIntersection[
       Ball[{Acen, Bcen, Ccen}, R1],
       Simplex[
        ElemPtsCoord[[  Delete[ClassifiedElements[[1, 1]], 1]  ]]
        ]
       ]
      ] != 0,

    ClassifiedElements[[1]],

    ClassifiedElements[[2]]
    ];
  ProjectedElementsIndex = 
   Union[Flatten[
     Take[ProjectedElementsIndex, Length[ProjectedElementsIndex], 1]]];

  ProjectedElementsIndex

  ]


ProjectedElemsIndex[ElemPtsCoord, ElemElemsMadeofPtsIndexandIdex, 
 62.84752406439229, 28.984094217780612, 72.55554474676572, 12.5, RetainElementsIndex]

Questions Background:

I have a cube sized 100*100*100, I have meshed it with fine meshes. So the mesh elements set is ElemElemsMadeofPtsIndex, The mesh points set is ElemPtsCoord.

And because ElemElemsMadeofPtsIndex and ElemPtsCoord do not tell us elements(points)' oder. So I add the points(elements) a order number to the first of the sublists:

ElemElemsMadeofPtsIndex-> ElemElemsMadeofPtsIndexandIdex

ElemPtsCoord-> ElemPtsCoordandIndex.

Math Problem:

A sphere(ball): the center is located at (Acen, Bcen, Ccen) and radius is R1).

and I want to find which elements belong(located in the ball or intersected with the ball) to the zone of the ball.

so the problem become to find sublists of list ElemElemsMadeofPtsIndexandIdex, which the volume of region which mesh tetrahedron and the ball intersected is not null.

$\endgroup$
  • $\begingroup$ possible duplicate of Find sublists of a list that match some conditions $\endgroup$ – RunnyKine Jul 22 '14 at 4:51
  • $\begingroup$ @ RunnyKine, No this is another questions. and I have not finished it yet, another question $\endgroup$ – YuYong Jul 22 '14 at 5:02
  • $\begingroup$ Okay, next time finish your questions before posting. It looked like a duplicate initially. $\endgroup$ – RunnyKine Jul 22 '14 at 5:05
  • $\begingroup$ @ RunnyKine Thanks,I am sorry $\endgroup$ – YuYong Jul 22 '14 at 5:07
0
$\begingroup$
    CountTrue[list_] := Module[
       {},
       Count[list, True] >= 1
       ];

    ProjectedElemsIndex[ElemPtsCoord_, ElemElemsMadeofPtsIndex_, Acen_, 
       Bcen_, Ccen_, R1_, RetainElementsIndex_] := Module[
       {RELPtsIndex, temp1, temp2},

       RELPtsIndex = ElemElemsMadeofPtsIndex[[RetainElementsIndex]];

       temp1 = 
        Reap[Do[Sow[  ElemPtsCoord[[RELPtsIndex[[i]]  ]] ], {i, 1, 
            Length[RELPtsIndex]} ] ][[2, 1]];
       temp2 = Table[Table[{Acen, Bcen, Ccen}, {4}], {Length[temp1]}];

       Flatten[
        Position[
         Map[
          CountTrue,
          Negative[
           Map[Norm, temp1 - temp2, {2}] - 
            Table[Table[R1, {4}], {Length[temp1 - temp2]}]
           ]
          ],
         True]
        ]

       ];

Length[ProjectedElemsIndex[ElemPtsCoord, ElemElemsMadeofPtsIndex, 40, 
   40, 40, 10, RetainElementsIndex]] // Timing
| improve this answer | |
$\endgroup$
  • $\begingroup$ This answer need to be optimistic $\endgroup$ – user18688 Jul 22 '14 at 12:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.