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.

Question2::

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]}
  ];

BoundaryPtsCoordandIndex = GatherBy[ElemPtsCoordandIndex,
   (#1[[2]] == 0) || (#1[[2]] == 100) || (#1[[3]] == 0) || (#1[[3]] ==
        100) || (#1[[4]] == 0) || (#1[[4]] == 100) &];

BoundaryPtsCoordandIndex = If[
   BoundaryPtsCoordandIndex[[1, 1]][[2]] == 0 ||
    BoundaryPtsCoordandIndex[[1, 1]][[2]] == 100 ||
    BoundaryPtsCoordandIndex[[1, 1]][[3]] == 0 ||
    BoundaryPtsCoordandIndex[[1, 1]][[3]] == 100 ||
    BoundaryPtsCoordandIndex[[1, 1]][[4]] == 0 ||
    BoundaryPtsCoordandIndex[[1, 1]][[4]] == 100,

   BoundaryPtsCoordandIndex[[1]],
   BoundaryPtsCoordandIndex[[2]]];


BoundaryPtsIndex = 
  Flatten[Take[BoundaryPtsCoordandIndex, 
    Length[BoundaryPtsCoordandIndex], 1]];

ClassifiedElements = 
  GatherBy[ElemElemsMadeofPtsIndexandIdex, 
   Intersection[BoundaryPtsIndex, Delete[#1, 1]] == {} &]

Introduction:

I have a cube sized 100*100*100. and I mesh it with fine mesh. So the total numbers of nodes and elements are huge.

and I have to find mesh tetrahedrons which are on the boundary.

First, Find points on the boundary faces-> "BoundaryPtsCoordandIndex".

Second, Find elements which has any point of "BoundaryPtsCoordandIndex".

This is my problem:

So the math problem is:

list1 = ElemElemsMadeofPtsIndexandIdex;

list2 = BoundaryPtsIndex;

and Find all sublists of list1 that contain any element of list2 ..?

share|improve this question

2 Answers 2

up vote 1 down vote accepted

This seems to work pretty well. Edit: now a bit faster.

rls = Thread[list2 -> True] // Dispatch;

Pick[list1, TrueQ /@ Or @@@ Replace[list1, rls, {2}]] // Length // Timing
{0.546003, 75230}

On my system this takes just over half a second to find 75,230 boundary sets.

(I am running version 10.0.0 under Windows.)


Looking at Simon's answer I find the whole killme thing extraneous. It is sufficient to merely replace all list2 elements with one that does not appear in anywhere in list1. I would instead write:

rls = Thread[list2 -> -1] // Dispatch;

Complement[list1, Replace[list1, rls, {2}]] // Length // Timing
{0.327602, 75230}

For me this also tests faster than Simon's code, which takes:

{0.374402, 75230}
share|improve this answer
    
@ Mr.Wizard, I am sorry to say what's the meaning of "/@ Or @@@", It seems that "TrueQ /@ Or @@@ Replace[list1, rls, {2}]" not equal to: "TrueQ /@ Replace[list1, rls, {2}]" Or "TrueQ /@ Replace[list1, rls, {2}]". I know that /@(Map) and @@@(Apply). –  YuYong Jul 22 at 1:32
    
@YuYong Or @@@ is to convert a list containing True to True: Or @@@ {{1, 2, 3}, {4, True, 6}} outputs {1 || 2 || 3, True}. Then TrueQ /@ converts anything not explicitly True` to False: TrueQ /@ {1 || 2 || 3, True} outputs {False, True}. –  Mr.Wizard Jul 22 at 2:12
    
@ Mr.Wizard,Yes, You are right. I think I have to Learn Mathematica's Functional Programming carefully. The speed of Conventional Programming is so slow. and I also have another similar questions: mathematica.stackexchange.com/questions/55533/… you help me again? –  YuYong Jul 22 at 5:28
    
Ah, of course! I feel silly now :-) –  Simon Woods Jul 22 at 12:10

On my system this is a little bit faster than Mr. Wizard's code, though it's a slightly odd way to approach the problem. I use a symbol which makes a list vanish when it appears as an element of that list...

killme /: {___, killme, ___} = Unevaluated[];

rls = Thread[list2 -> killme] // Dispatch;

result = Complement[list1, Replace[list1, rls, {2}]];
share|improve this answer
    
Yes, It seems so strange, because that the results of Replace[list1, rls, {2}] seems the same as list1 ... But the idea is great. and easy to understand. Cool! –  YuYong Jul 22 at 1:28
    
Confirmed that this is a bit faster. Let me see if I can beat it! –  Mr.Wizard Jul 22 at 2:15

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.