← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/mesh/GTSInterface.cpp

  • Committer: Johannes Ring
  • Date: 2008-03-05 22:43:06 UTC
  • Revision ID: johannr@simula.no-20080305224306-2npsdyhfdpl2esji
The BIG commit!

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/kernel/mesh/GTSInterface.cpp
1
 
// Copyright (C) 2006 Anders Logg.
2
 
// Licensed under the GNU LGPL Version 2.1.
3
 
//
4
 
// Modified by Johan Jansson 2006.
5
 
// Modified by Ola Skavhaug 2006.
6
 
// Modified by Dag Lindbo 2008.
7
 
//
8
 
// First added:  2006-06-21
9
 
// Last changed: 2006-12-01
10
 
 
11
 
#include <dolfin/dolfin_log.h>
12
 
#include <dolfin/Array.h>
13
 
#include <dolfin/Mesh.h>
14
 
#include <dolfin/Facet.h>
15
 
#include <dolfin/Vertex.h>
16
 
#include <dolfin/Cell.h>
17
 
 
18
 
#include <dolfin/GTSInterface.h>
19
 
 
20
 
#ifdef HAVE_GTS_H
21
 
#include <gts.h>
22
 
#endif
23
 
 
24
 
using namespace dolfin;
25
 
 
26
 
//-----------------------------------------------------------------------------
27
 
GTSInterface::GTSInterface(Mesh& m) : mesh(m), tree(0) 
28
 
{
29
 
  buildCellTree();
30
 
}
31
 
//-----------------------------------------------------------------------------
32
 
GTSInterface::~GTSInterface() 
33
 
{
34
 
#ifdef HAVE_GTS_H
35
 
  gts_bb_tree_destroy(tree, 1);
36
 
#endif
37
 
}
38
 
//-----------------------------------------------------------------------------
39
 
GtsBBox* GTSInterface::bboxCell(Cell& c)
40
 
{
41
 
#ifdef HAVE_GTS_H
42
 
  GtsBBox* bbox;
43
 
  Point p;
44
 
 
45
 
  VertexIterator v(c);
46
 
  p = v->point();
47
 
 
48
 
  bbox = gts_bbox_new(gts_bbox_class(), (void *)c.index(),
49
 
                      p.x(), p.y(), p.z(),
50
 
                      p.x(), p.y(), p.z());
51
 
  
52
 
  for(++v; !v.end(); ++v)
53
 
    {
54
 
      p = v->point();
55
 
      if (p.x() > bbox->x2) bbox->x2 = p.x();
56
 
      if (p.x() < bbox->x1) bbox->x1 = p.x();
57
 
      if (p.y() > bbox->y2) bbox->y2 = p.y();
58
 
      if (p.y() < bbox->y1) bbox->y1 = p.y();
59
 
      if (p.z() > bbox->z2) bbox->z2 = p.z();
60
 
      if (p.z() < bbox->z1) bbox->z1 = p.z();
61
 
    }
62
 
  return bbox;
63
 
 
64
 
#else
65
 
  error("missing GTS");
66
 
  return 0;
67
 
#endif
68
 
}
69
 
//-----------------------------------------------------------------------------
70
 
GtsBBox* GTSInterface::bboxPoint(const Point& p)
71
 
{
72
 
#ifdef HAVE_GTS_H
73
 
 
74
 
  GtsBBox* bbox;
75
 
 
76
 
  bbox = gts_bbox_new(gts_bbox_class(), (void *)0,
77
 
                      p.x(), p.y(), p.z(),
78
 
                      p.x(), p.y(), p.z());
79
 
  
80
 
  return bbox;
81
 
 
82
 
#else
83
 
  error("missing GTS");
84
 
  return 0;
85
 
#endif
86
 
}
87
 
//-----------------------------------------------------------------------------
88
 
GtsBBox* GTSInterface::bboxPoint(const Point& p1, const Point& p2)
89
 
{
90
 
#ifdef HAVE_GTS_H
91
 
 
92
 
  GtsBBox* bbox;
93
 
 
94
 
  real x1, x2; 
95
 
  real y1, y2;
96
 
  real z1, z2;
97
 
 
98
 
  if(p1.x()<p2.x()){
99
 
    x1 = p1.x();
100
 
    x2 = p2.x(); }
101
 
  else {
102
 
    x1 = p2.x();
103
 
    x2 = p1.x(); }
104
 
  if(p1.y()<p2.y()){
105
 
    y1 = p1.y();
106
 
    y2 = p2.y(); }
107
 
  else {
108
 
    y1 = p2.y();
109
 
    y2 = p1.y(); }
110
 
  if(p1.z()<p2.z()){
111
 
    z1 = p1.z();
112
 
    z2 = p2.z(); }
113
 
  else {
114
 
    z1 = p2.z();
115
 
    z2 = p1.z(); }
116
 
 
117
 
  bbox = gts_bbox_new(gts_bbox_class(), (void *)0,
118
 
                      x1, y1, z1,
119
 
                      x2, y2, z2);
120
 
  
121
 
  return bbox;
122
 
 
123
 
#else
124
 
  error("missing GTS");
125
 
  return 0;
126
 
#endif
127
 
}
128
 
//-----------------------------------------------------------------------------
129
 
void GTSInterface::buildCellTree()
130
 
{
131
 
#ifdef HAVE_GTS_H
132
 
 
133
 
  if(tree)
134
 
    warning("tree already initialized");
135
 
 
136
 
  GSList* bboxes = NULL;
137
 
 
138
 
  for(CellIterator ci(mesh); !ci.end(); ++ci)
139
 
    {
140
 
      Cell& c = *ci;
141
 
      bboxes = g_slist_prepend(bboxes, bboxCell(c));
142
 
    }
143
 
  
144
 
  tree = gts_bb_tree_new(bboxes);
145
 
  g_slist_free(bboxes);
146
 
 
147
 
#else
148
 
  error("missing GTS");
149
 
#endif
150
 
}
151
 
//-----------------------------------------------------------------------------
152
 
void GTSInterface::overlap(Cell& c, Array<uint>& cells)
153
 
{
154
 
#ifdef HAVE_GTS_H
155
 
  GtsBBox* bbprobe;
156
 
  GtsBBox* bb;
157
 
  GSList* overlaps = 0, *overlaps_base;
158
 
  uint boundedcell;
159
 
 
160
 
  CellType& type = mesh.type();
161
 
 
162
 
  bbprobe = bboxCell(c);
163
 
 
164
 
  overlaps = gts_bb_tree_overlap(tree, bbprobe);
165
 
  overlaps_base = overlaps;
166
 
 
167
 
  while(overlaps)
168
 
    {
169
 
      bb = (GtsBBox *)overlaps->data;
170
 
      boundedcell = (uint)(long)bb->bounded;
171
 
 
172
 
      Cell close(mesh, boundedcell);
173
 
 
174
 
      if(type.intersects(c, close))
175
 
        {
176
 
          cells.push_back(boundedcell);
177
 
        }
178
 
      overlaps = overlaps->next;
179
 
    }
180
 
  
181
 
  g_slist_free(overlaps_base);
182
 
  gts_object_destroy(GTS_OBJECT(bbprobe));
183
 
 
184
 
#else
185
 
  error("missing GTS");
186
 
#endif
187
 
}
188
 
//-----------------------------------------------------------------------------
189
 
void GTSInterface::overlap(Point& p, Array<uint>& cells)
190
 
{
191
 
#ifdef HAVE_GTS_H
192
 
  GtsBBox* bbprobe;
193
 
  GtsBBox* bb;
194
 
  GSList* overlaps = 0, *overlaps_base;
195
 
  uint boundedcell;
196
 
 
197
 
  CellType& type = mesh.type();
198
 
 
199
 
  bbprobe = bboxPoint(p);
200
 
 
201
 
  overlaps = gts_bb_tree_overlap(tree, bbprobe);
202
 
  overlaps_base = overlaps;
203
 
 
204
 
  while(overlaps)
205
 
    {
206
 
      bb = (GtsBBox *)overlaps->data;
207
 
      boundedcell = (uint)(long)bb->bounded;
208
 
 
209
 
      Cell close(mesh, boundedcell);
210
 
 
211
 
      if(type.intersects(close, p))
212
 
        cells.push_back(boundedcell);
213
 
 
214
 
      overlaps = overlaps->next;
215
 
    }
216
 
 
217
 
  g_slist_free(overlaps_base);
218
 
  gts_object_destroy(GTS_OBJECT(bbprobe));
219
 
 
220
 
#else
221
 
  error("missing GTS");
222
 
#endif
223
 
}
224
 
//-----------------------------------------------------------------------------
225
 
void GTSInterface::overlap(Point& p1, Point& p2, Array<uint>& cells)
226
 
{
227
 
#ifdef HAVE_GTS_H
228
 
  GtsBBox* bbprobe;
229
 
  GtsBBox* bb;
230
 
  GSList* overlaps = 0,*overlaps_base;
231
 
  uint boundedcell;
232
 
 
233
 
  bbprobe = bboxPoint(p1,p2);
234
 
 
235
 
  overlaps = gts_bb_tree_overlap(tree, bbprobe);
236
 
  overlaps_base = overlaps;
237
 
 
238
 
  while(overlaps)
239
 
    {
240
 
      bb = (GtsBBox *)overlaps->data;
241
 
      boundedcell = (uint)(long)bb->bounded;
242
 
 
243
 
      cells.push_back(boundedcell);
244
 
 
245
 
      overlaps = overlaps->next;
246
 
    }
247
 
  g_slist_free(overlaps_base);
248
 
  gts_object_destroy(GTS_OBJECT(bbprobe));
249
 
 
250
 
#else
251
 
  error("missing GTS");
252
 
#endif
253
 
}
254
 
//-----------------------------------------------------------------------------