Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.
  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.

share|improve this question
    
possible duplicate of Find sublists of a list that match some conditions –  RunnyKine Jul 22 at 4:51
    
@ RunnyKine, No this is another questions. and I have not finished it yet, another question –  YuYong Jul 22 at 5:02
    
Okay, next time finish your questions before posting. It looked like a duplicate initially. –  RunnyKine Jul 22 at 5:05
    
@ RunnyKine Thanks,I am sorry –  YuYong Jul 22 at 5:07

1 Answer 1

up vote 0 down vote accepted
    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
share|improve this answer
    
This answer need to be optimistic –  user18688 Jul 22 at 12:52

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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