# source:branches/parallel/Ipopt/contrib/PetscInterface/ProblemGeometry.C@1776

Last change on this file since 1776 was 1776, checked in by andreasw, 3 years ago

experiments with AMR

• Property svn:eol-style set to `native`
• Property svn:keywords set to `Author Date Id Revision`
File size: 48.1 KB
Line
3// This code is published under the Common Public License.
4//
5// \$Id\$
6//
7// Authors:  Johannes Huber     IBM        2010-09-03
8
9#include "ProblemGeometry.h"
10#include <stdexcept>
11
12///////////////////////////////////////////////////////////////////////////////////
13//  ProblemGeometry::Item
14ProblemGeometry::Item::Item(const std::vector<double>& min, const std::vector<double>& max)
15{
16  assert(min.size()==max.size());
17  min_=min;
18  max_=max;
19  ValidateMinMax();
20  BoundaryMarker.resize(2*min_.size(),-1); // those are for floor and ceiling
21}
22
23void ProblemGeometry::Item::ValidateMinMax()
24{
25  for (int i=0;i<min_.size();i++) {
26    if (min_[i]>max_[i]) {
27      double tmp=min_[i];
28      min_[i]=max_[i];
29      max_[i]=tmp;
30    }
31  }
32}
33
34bool ProblemGeometry::Item::IsInArea(const std::vector<double>& pt, double epsilon/*=1e-8*/)
35{
36  assert(pt.size()==min_.size());
37  int dim = min_.size();
38  for ( int i=0; i<dim; i++ ) {
39    if (min_[i]-epsilon>pt[i])
40      return false;
41    if (max_[i]+epsilon<pt[i])
42      return false;
43  }
44  return true;
45}
46
47int ProblemGeometry::Item::GetClosestBoundary(const std::vector<double>& pt)
48{
49  assert(pt.size()==min_.size());
50  int dim = min_.size();
51  double min_dist=fabs(pt[0]-min_[0]);
52  int cur_min_dir=0;
53  for (int i=0; i<dim; i++) {
54    if (fabs(pt[i]-min_[i])<min_dist) {
55      min_dist = fabs(pt[i]-min_[i]);
56      cur_min_dir = 2*i;
57    }
58    if (fabs(max_[i]-pt[i])<min_dist) {
59      min_dist = fabs(max_[i]-pt[i]);
60      cur_min_dir = 2*i+1;
61    }
62  }
63  return cur_min_dir;
64}
65
66int ProblemGeometry::Item::GetBoundaryMarker(const std::vector<double>& pt)
67{
68  assert(pt.size()==min_.size());
69  if (!IsInArea(pt)) return -1;
70  return BoundaryMarker[GetClosestBoundary(pt)];
71}
72
73//////////////////////////////////////////////////////////////////////////////////////
74// ProblemGeometry
75ProblemGeometry::ProblemGeometry()
76{
77  BoundaryCondition BC(0.0,1.0,0.0,0.0,1.0,0.0);  // 2 * homogene Neumann
78  _BoundCond.push_back(BC); //   Add one Dummy to start from 1 (0:Interior Point)
79  _BoundCond.push_back(BC);
80  NextFreeBoundaryMarker = 2;
81  _h = 0.0;
82}
83
84void ProblemGeometry::AddEquipment(const std::vector<double>& min, const std::vector<double>& max, double TempEquip, double kappa)
85{
86  assert(min.size()==GetDim());
87  assert(max.size()==GetDim());
88  Item item(min,max);
89  assert(_BoundCond.size()==NextFreeBoundaryMarker);
90  int BoundMark = NextFreeBoundaryMarker;
91  // 1.st & 3rd Boundary: Head exchange
92  item.BoundaryMarker[2] = BoundMark;
93  item.BoundaryMarker[3] = BoundMark;
94  // Phi: hom. Neum. T: dT/dn = Kappa(T-Teq) => 1*dT/dn = Kappa*T -Kappa*Teq
95  _BoundCond.push_back(BoundaryCondition(0.0,1.0,0.0,kappa,-1.0,-kappa*TempEquip));
96  // No control parameters
97  // 2.nd Boundary: Isolation
98  item.BoundaryMarker[0] = BoundMark+1;
99  item.BoundaryMarker[1] = BoundMark+1;
100  // Phi: hom. Neum. T: hom. Neum
101  _BoundCond.push_back(BoundaryCondition(0.0,1.0,0.0,0.0,1.0,0.0));
102  // No control parameters
103  NextFreeBoundaryMarker+=2;
104  if (GetDim()>2) { // Equipmentceiling: heat exchange (same as on side)
105    item.BoundaryMarker[4] = -1;   // should never be accessed, since Equipment is standing on floor, so no boundary
106    int BoundMark = NextFreeBoundaryMarker++;
107    // 1.st & 3rd Boundary: Head exchange
108    item.BoundaryMarker[5] = BoundMark;
109    // Phi: hom. Neum. T: dT/dn = Kappa(T-Teq)
110    _BoundCond.push_back(BoundaryCondition(0.0,1.0,0.0,kappa,1.0,kappa*TempEquip));
111  }
112  _Equip.push_back(item);
113}
114
115void ProblemGeometry::AddAC(const std::vector<double>& min, const std::vector<double>& max, double vAc, double TempAc)
116{
117  assert(min.size()==GetDim());
118  assert(max.size()==GetDim());
119  _ItemAtWall.clear();  // invalidate item at wall and rebuild it, when the mesh is created
120  Item item(min,max);
121  int BoundMark = NextFreeBoundaryMarker++;
122  for (unsigned int iBd=0;iBd<2*GetDim();iBd++)
123    item.BoundaryMarker[iBd] = BoundMark;
124  assert(_BoundCond.size()==BoundMark);
125  _BoundCond.push_back(BoundaryCondition(0.0,-1.0,vAc,1.0,0.0,-TempAc)); // dPhi/dn = -vAc,T = TAc, i.e. 0 dT/dn = 1*T-TAc
126  _ParamIdx2BCParam.push_back(ControlParameter(BoundMark,2));
127  //_ParamIdx2BCParam.push_back(ControlParameter(BoundMark,5));
128  _AC.push_back(item);
129}
130
131#ifndef EXHAUST_AS_CONTROL
132void ProblemGeometry::AddExhaust(const std::vector<double>& min, const std::vector<double>& max)
133{
134  assert(min.size()==GetDim());
135  assert(max.size()==GetDim());
136  _ItemAtWall.clear();
137  Item item(min,max);
138  unsigned int BoundMark = NextFreeBoundaryMarker++;
139  for ( int iBd=0;iBd<2*GetDim();iBd++)
140    item.BoundaryMarker[iBd] = BoundMark;
141  assert(_BoundCond.size()==BoundMark);
142  _BoundCond.push_back(BoundaryCondition(1.0,0.0,0.0,0.0,1.0,0.0));         // Phi = 0, dT/dn = 0
143  // No Control Parameter
144  _Exh.push_back(item);
145}
146#else
147void ProblemGeometry::AddExhaust(const std::vector<double>& min, const std::vector<double>& max)
148{
149  assert(min.size()==GetDim());
150  assert(max.size()==GetDim());
151  _ItemAtWall.clear();
152  Item item(min,max);
153  unsigned int BoundMark = NextFreeBoundaryMarker++;
154  for ( int iBd=0;iBd<2*GetDim();iBd++)
155    item.BoundaryMarker[iBd] = BoundMark;
156  assert(_BoundCond.size()==BoundMark);
157  // TODO: Retried vExh from input
158  double vExh = 1.0;
159  double TempExh = 1e30;
160  _BoundCond.push_back(BoundaryCondition(0.0,+1.0,vExh,1.0,0.0,-TempExh)); // dPhi/dn = -vAc,T = TAc, i.e. 0 dT/dn = 1*T-TAc
161  _ParamIdx2BCParam.push_back(ControlParameter(BoundMark,2));
162  //_ParamIdx2BCParam.push_back(ControlParameter(BoundMark,5));
163  _Exh.push_back(item);
164}
165#endif
166
168{
169  static char Buf[1024];
170  int n;
171  double Vals[8];
172  _h=0;
173  while (!is.eof()) {
174    is >> Buf;
175    if (strlen(Buf)<1)
176      continue;
177    if (Buf[0]=='#')
178      continue;
179    n = sscanf(Buf,"h=%lf",Vals);
180    if ( (n==1) ) {
181      _h = Vals[0];
182      break;
183    }
184    else {
185      std::cout << "WARNING while reading input file and looking for \"h=*\": Line in wrong position?" << std::endl;
186      std::cout << Buf << std::endl;
187      continue;
188    }
189  }
190
191  _RoomSize.clear();
192  while (!is.eof()) {
193    is >> Buf;
194    if (strlen(Buf)<1)
195      continue;
196    if (Buf[0]=='#')
197      continue;
198    n = sscanf(Buf,"Roomsize=%lf,%lf,%lf",Vals,Vals+1,Vals+2);
199    if ( (n==2) || (n==3) ) {
200      for (int i=0;i<n;i++)
201        _RoomSize.push_back(Vals[i]);
202      break;
203    }
204    else {
205      std::cout << "WARNING while reading input file and looking for \"Roomsize=*,*[,*]\": Line in wrong position?" << std::endl;
206      std::cout << Buf << std::endl;
207      continue;
208    }
209  }
210  if (0==_RoomSize.size()) {
211    throw std::runtime_error("Can't find line \"Roomsize=*,*[,*]\"");
212  }
213  else {
214    std::cout << "Roomsize=" << _RoomSize[0] <<"," << _RoomSize[1];
215    if (_RoomSize.size()>2)
216      std::cout << "," << _RoomSize[2];
217    std::cout << std::endl;
218  }
219
220  _Equip.clear();
221  _AC.clear();
222  _Exh.clear();
223  std::vector<double> tmp_min;
224  std::vector<double> tmp_max;
225  while (!is.eof()) {
226    tmp_max.clear();
227    tmp_min.clear();
228    is >> Buf;
229    std::cout << Buf << "line read" << std::endl;
230    if (strlen(Buf)<1)
231      continue;
232    if (Buf[0]=='#')
233      continue;
234
235    n = sscanf(Buf,"AC:Min=%lf,%lf,%lf;",Vals,Vals+1,Vals+2);
236    if (n>1) {
237      if (n==2) {
238        n = sscanf(Buf,"AC:Min=%lf,%lf;Max=%lf,%lf;vAC=%lf;TAC=%lf",Vals,Vals+1,Vals+2,Vals+3,Vals+4,Vals+5);
239        if (n!=6) {
240          std::string str("Can't interpret line as AC in 2D:");
241          str += Buf;
242          throw std::runtime_error(str);
243        }
244        tmp_min.push_back(Vals[0]);
245        tmp_min.push_back(Vals[1]);
246        tmp_max.push_back(Vals[2]);
247        tmp_max.push_back(Vals[3]);
249      }
250      else {
251        n = sscanf(Buf,"AC:Min=%lf,%lf,%lf;Max=%lf,%lf,%lf;vAC=%lf;TAC=%lf",Vals,Vals+1,Vals+2,Vals+3,Vals+4,Vals+5,Vals+6,Vals+7);
252        if (n!=8) {
253          std::string str("Can't interpret line as AC in 3D:");
254          str += Buf;
255          throw std::runtime_error(str);
256        }
257        tmp_min.push_back(Vals[0]);
258        tmp_min.push_back(Vals[1]);
259        tmp_min.push_back(Vals[2]);
260        tmp_max.push_back(Vals[3]);
261        tmp_max.push_back(Vals[4]);
262        tmp_max.push_back(Vals[5]);
264      }
265    }
266    else {
267      n = sscanf(Buf,"Exh:Min=%lf,%lf,%lf;",Vals,Vals+1,Vals+2);
268      if (n>1) {
269        if (n==2) {
270          n = sscanf(Buf,"Exh:Min=%lf,%lf;Max=%lf,%lf",Vals,Vals+1,Vals+2,Vals+3);
271          if (n!=4) {
272            std::string str("Can't interpret line as Exh in 2D:");
273            str += Buf;
274            throw std::runtime_error(str);
275          }
276          tmp_min.push_back(Vals[0]);
277          tmp_min.push_back(Vals[1]);
278          tmp_max.push_back(Vals[2]);
279          tmp_max.push_back(Vals[3]);
281        }
282        else {
283          n = sscanf(Buf,"Exh:Min=%lf,%lf,%lf;Max=%lf,%lf,%lf",Vals,Vals+1,Vals+2,Vals+3,Vals+4,Vals+5);
284          if (n!=6) {
285            std::string str("Can't interpret line as Exh in 3D:");
286            str += Buf;
287            throw std::runtime_error(str);
288          }
289          tmp_min.push_back(Vals[0]);
290          tmp_min.push_back(Vals[1]);
291          tmp_min.push_back(Vals[2]);
292          tmp_max.push_back(Vals[3]);
293          tmp_max.push_back(Vals[4]);
294          tmp_max.push_back(Vals[5]);
296        }
297      }
298      else {
299        n = sscanf(Buf,"Equip:Min=%lf,%lf,%lf;",Vals,Vals+1,Vals+2);
300        if (n>1) {
301          if (n==2) {
302            n = sscanf(Buf,"Equip:Min=%lf,%lf;Max=%lf,%lf;TEquip=%lf;kappa=%lf",Vals,Vals+1,Vals+2,Vals+3,Vals+4,Vals+5);
303            if (n!=6) {
304              std::string str("Can't interpret line as Equip in 2D:");
305              str += Buf;
306              throw std::runtime_error(str);
307            }
308            tmp_min.push_back(Vals[0]);
309            tmp_min.push_back(Vals[1]);
310            tmp_max.push_back(Vals[2]);
311            tmp_max.push_back(Vals[3]);
313          }
314          else {
315            n = sscanf(Buf,"Equip:Min=%lf,%lf,%lf;Max=%lf,%lf,%lf;TEquip=%lf;kappa=%lf",Vals,Vals+1,Vals+2,Vals+3,Vals+4,Vals+5,Vals+6,Vals+7);
316            if (n!=8) {
317              std::string str("Can't interpret line as Equip in 3D:");
318              str += Buf;
319              throw std::runtime_error(str);
320            }
321            tmp_min.push_back(Vals[0]);
322            tmp_min.push_back(Vals[1]);
323            tmp_min.push_back(Vals[2]);
324            tmp_max.push_back(Vals[3]);
325            tmp_max.push_back(Vals[4]);
326            tmp_max.push_back(Vals[5]);
328          }
329        }
330        else {
331          std::string str("Can't interpret line:");
332          str += Buf;
333          throw std::runtime_error(str);
334        }
335      }
336    }
337  }
338  std::cout << "Input file read" << std::endl;
339}
340
341int ProblemGeometry::GetBoundaryMarker(const std::vector<double>& pt)
342{
343  int iWall = GetWall(pt);
344  int Tmp;
345  if (iWall<0)  {
346    for (int iEq=0;iEq<_Equip.size();iEq++) {
347      Tmp=_Equip[iEq].GetBoundaryMarker(pt);
348      if (Tmp>=0)
349        return Tmp;
350    }
351    return -1;
352  }
353  else {
354    if (iWall>3) // floor and ceiling
355      return 1;
356    ValidateItemAtWall();
357    std::list< Item* >::iterator ppItem;
358    ppItem = _ItemAtWall[iWall].begin();
359    for (;ppItem!=_ItemAtWall[iWall].end();ppItem++) {
360      Tmp=(*ppItem)->GetBoundaryMarker(pt);
361      if (Tmp>=0)
362        return Tmp;
363    }
364    return 1;
365  }
366}
367
368int ProblemGeometry::GetWall(const std::vector<double>& pt)
369{
370  assert(pt.size()==GetDim());
371  std::vector<double> Null;
372  Null.resize(GetDim(),0.0);
373
374  Item room(Null,_RoomSize);
375  if (room.IsInArea(pt,-1e-8)) // Check if it's inside
376    return -1;
377  if (room.IsInArea(pt))
378    return room.GetClosestBoundary(pt);
379  else
380    return -1;
381}
382
383// build up _TimeAtWall member for each 4 walls
384void ProblemGeometry::ValidateItemAtWall()
385{
386  // _ItemAtWall: Index: 0:x0=0, 1:x0=Max, 2:x1=0, 3:x1=Max
387  double epsilon = 1e-8;
388  if (_ItemAtWall.size()>0)
389    return;
390  std::list< Item* > EmptyList;
391  _ItemAtWall.resize(4,EmptyList);
392  Item* pItem;
393  unsigned int iItem;
394  for (iItem=0;iItem<_AC.size();iItem++) {
395    pItem = &(_AC[iItem]);
396    for (int iDim=0;iDim<2;iDim++) // No check for iDim=2, thats floor or ceiling, special case
397    {
398      if ( fabs(pItem->min_[iDim])<epsilon && fabs(pItem->max_[iDim])<epsilon )
399      {
400        _ItemAtWall[2*iDim].push_back(pItem);
401        continue;
402      }
403      if ( fabs(pItem->min_[iDim]-_RoomSize[iDim])<epsilon && fabs(pItem->max_[iDim]-_RoomSize[iDim])<epsilon )
404      {
405        _ItemAtWall[2*iDim+1].push_back(pItem);
406        continue;
407      }
408    }
409  }
410
411  for (iItem=0;iItem<_Exh.size();iItem++) {
412    pItem = &(_Exh[iItem]);
413    for (int iDim=0;iDim<2;iDim++) // No check for iDim=2, thats floor or ceiling, special case
414    {
415      if ( fabs(pItem->min_[iDim])<epsilon && fabs(pItem->max_[iDim])<epsilon )
416      {
417        _ItemAtWall[2*iDim].push_back(pItem);
418        continue;
419      }
420      if ( fabs(pItem->min_[iDim]-_RoomSize[iDim])<epsilon && fabs(pItem->max_[iDim]-_RoomSize[iDim])<epsilon )
421      {
422        _ItemAtWall[2*iDim+1].push_back(pItem);
423        continue;
424      }
425    }
426  }
427}
428
429
430// return index of wall at which item is located
431int ProblemGeometry::GetWallItemWall(const Item& item)
432{
433  // Returnvalue: Min x1 -> 0, Max x1->1, Min x2 -> 2 Max x2 -> 3, Min x3 -> 4 Max x3 -> 5
434  double epsilon = 1e-8;
435  for (int iDim=0;iDim<GetDim();iDim++) {
436    if ( (fabs(item.min_[iDim])<epsilon) && (fabs(item.max_[iDim])<epsilon) )
437      return 2*iDim;
438    if ( (fabs(item.min_[iDim]-_RoomSize[iDim])<epsilon) && (fabs(item.max_[iDim]-_RoomSize[iDim])<epsilon) )
439      return 2*iDim+1;
440  }
441  // We should never get here
442  assert(false);
443  return -1;
444}
445
446// Update p_mesh to include BoundaryMarkers for all nodes that have Dirichlet
447void ProblemGeometry::SetBoundaryInfo(libMesh::Mesh* p_mesh)
448{
449  const double eps = 1e-8;
450//  p_mesh->boundary_info->clear();
451
452  //mark element sides
453  std::vector<double> Pos;
454  Pos.resize(GetDim());
455  MeshBase::const_element_iterator el = p_mesh->active_local_elements_begin();
456  const MeshBase::const_element_iterator end_el = p_mesh->active_local_elements_end();
457  // loop over all active local elements
458  for ( ; el != end_el; ++el) {
459    for (unsigned int side=0; side<(*el)->n_sides(); side++) {
460      if ((*el)->neighbor(side) == NULL) {
461        // this is element side on boundary
462        libMesh::Point Center;
463        for (unsigned int iPt=0; iPt<(*el)->n_nodes();iPt++)
464          if ((*el)->is_node_on_side(iPt,side))
465            Center += (*el)->point(iPt);
466        Center = (1.0/GetDim()) * Center;
467        for (unsigned int iDim=0;iDim<GetDim();iDim++)
468          Pos[iDim] = Center(iDim);
469        int BoundaryMarker = GetBoundaryMarker(Pos);
470        assert(BoundaryMarker!=-1);
472
473        // mark nodes if Dirichlet cond.: Those are set at points (not sides) by putting diagonals into the matrix
474        if ( (fabs(_BoundCond[BoundaryMarker].PhiNeumannCoef)<eps) || (fabs(_BoundCond[BoundaryMarker].TNeumannCoef)<eps) ) {
475          // mark side points
476          for (unsigned int iPt=0; iPt<(*el)->n_nodes();iPt++)
477            if ((*el)->is_node_on_side(iPt,side))
479        }
480      }
481    }
482  }
483}
484
485/////////////////////////////////////////////////////////////////////
486//  3D
487void ProblemGeometry::SetFacette(tetgenio::facet *f, int Base, int Offset0, int Offset1, int Offset2, int Offset3, int nHoles)
488{
489  tetgenio::init(f);
490  f->numberofholes=nHoles;
491  if (nHoles>0)
492    f->holelist = new REAL[f->numberofholes * 3];
493  f->numberofpolygons=1+nHoles;
494  f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
495  tetgenio::polygon *p = f->polygonlist;
496  tetgenio::init(p);
497  p->numberofvertices = 4;
498  p->vertexlist = new int[p->numberofvertices];
499  p->vertexlist[0]=Base+Offset0;
500  p->vertexlist[1]=Base+Offset1;
501  p->vertexlist[2]=Base+Offset2;
502  p->vertexlist[3]=Base+Offset3;
503}
504
505void ProblemGeometry::SetPoint(REAL* pointlist,int idx,REAL x,REAL y,REAL z)
506{
507  pointlist[3*idx+0]=x;
508  pointlist[3*idx+1]=y;
509  pointlist[3*idx+2]=z;
510}
511
512
513void ProblemGeometry::Tetgen2Mesh(const tetgenio& tet, libMesh::UnstructuredMesh* p_mesh)
514{
515  const double eps(1e-8);
516  assert(GetDim()==3);
517  int dim = GetDim();
518  p_mesh->set_mesh_dimension(3);
520  unsigned long int iPt;
521  int BoundaryMarker;
522  for (iPt=0;iPt<tet.numberofpoints;iPt++)
524  // elements
525  libMesh::Elem *pElem;
526  for (unsigned int iEl=0; iEl<tet.numberoftetrahedra; iEl++) {
527    pElem = new libMesh::Tet4;
528    for (iPt=0;iPt<4;iPt++) {
529      pElem->set_node(iPt) = p_mesh->node_ptr(tet.tetrahedronlist[4*iEl+iPt]);
530    }
531    // Finally, add this element to the mesh
533  }
534
535  if (tet.neighborlist) {
536    for (unsigned int iEl=0; iEl<tet.numberoftetrahedra; iEl++) {
537      for (int iNeig=0;iNeig<dim+1;++iNeig) {
538        if (tet.neighborlist[4*iEl+iNeig]==-1) {
539          libMesh::Tet4 *pElem = static_cast<libMesh::Tet4*>(p_mesh->elem(iEl));
540          int iBdNode1 = tet.tetrahedronlist[(dim+1)*iEl+((iNeig+1)%(dim+1))];
541          int iBdNode2 = tet.tetrahedronlist[(dim+1)*iEl+((iNeig+2)%(dim+1))];
542          int iBdNode3 = tet.tetrahedronlist[(dim+1)*iEl+((iNeig+3)%(dim+1))];
543          int sum = iBdNode1+iBdNode2+iBdNode3;
544          for (int iSide=0;iSide<dim+1;iSide++) {
545            // sum is correct -> sid is correct (holds only, if side has one node less than elem)
546            if ( (pElem->node(pElem->side_nodes_map[iSide][0])
547                  +pElem->node(pElem->side_nodes_map[iSide][1])
548                  +pElem->node(pElem->side_nodes_map[iSide][2])) == sum ) {
549              std::vector<double> Center;
550              Center.resize(dim);
551              for (int iDim=0;iDim<dim;iDim++)
552                Center[iDim] = (tet.pointlist[dim*iBdNode1+iDim] + tet.pointlist[dim*iBdNode2+iDim]+ tet.pointlist[dim*iBdNode3+iDim])/dim;
553              int BoundaryMarker = GetBoundaryMarker(Center);
554              assert(BoundaryMarker!=-1);
556              if ( (fabs(_BoundCond[BoundaryMarker].PhiNeumannCoef)<eps) ||
557                   (fabs(_BoundCond[BoundaryMarker].TNeumannCoef)<eps) ) {
561              }
562            }
563          }
564        }
565      }
566    }
567  }
568}
569
570void ProblemGeometry::CreateMesh(libMesh::UnstructuredMesh* p_mesh)
571{
572  switch (GetDim()) {
573  case 2:
574    CreateMesh2D(p_mesh);
575    break;
576  case 3:
577    CreateMesh3D(p_mesh);
578    break;
579  default:
580    std::cout << "Dimension not implemented" << std::endl;
581    exit(0);
582  }
583}
584
585void ProblemGeometry::CreateMesh3D(libMesh::UnstructuredMesh* p_mesh)
586{
587  tetgenio tetgen_in;
588  //tetgenmesh      tetgen_mesh;
589
590  tetgen_in.mesh_dim = 3;  // Three-dimemsional accoordinates.
591  tetgen_in.numberofpointattributes = 0;  // no point attribute.
592  tetgen_in.firstnumber = 0;
593  tetgen_in.numberofcorners = 1;
594  tetgenio::facet *FacWall;
595  tetgenio::polygon *p;
596
597  tetgen_in.numberofpoints = 8+8*_Equip.size()+4*(_AC.size()+_Exh.size());
598  tetgen_in.pointlist = new REAL[3*tetgen_in.numberofpoints];
599  int iPtInList = 0;
600  int iFacInList = 0;
601
602  SetPoint(tetgen_in.pointlist,iPtInList++,0             ,0             ,0             );
603  SetPoint(tetgen_in.pointlist,iPtInList++,0             ,0             ,_RoomSize[2]);
604  SetPoint(tetgen_in.pointlist,iPtInList++,0             ,_RoomSize[1],0             );
605  SetPoint(tetgen_in.pointlist,iPtInList++,0             ,_RoomSize[1],_RoomSize[2]);
606  SetPoint(tetgen_in.pointlist,iPtInList++,_RoomSize[0],0             ,0             );
607  SetPoint(tetgen_in.pointlist,iPtInList++,_RoomSize[0],0             ,_RoomSize[2]);
608  SetPoint(tetgen_in.pointlist,iPtInList++,_RoomSize[0],_RoomSize[1],0             );
609  SetPoint(tetgen_in.pointlist,iPtInList++,_RoomSize[0],_RoomSize[1],_RoomSize[2]);
610
611  // Where are the AC's and the Exhausts?
612  ValidateItemAtWall();
613  std::list<Item*>::iterator ppItem;
614
615  int NumOfWallItems[4] = {0,0,0,0};
616  unsigned int iItem;
617  for (iItem=0;iItem<_AC.size();iItem++)
618    NumOfWallItems[GetWallItemWall(_AC[iItem])]++;
619
620  for (iItem=0;iItem<_Exh.size();iItem++)
621    NumOfWallItems[GetWallItemWall(_Exh[iItem])]++;
622
623  // facet marker:0: normal boundary condition, rest increasing
624  int iFacMakerUsed = 1;
625  tetgen_in.numberoffacets = 6+5*_Equip.size()+_AC.size()+_Exh.size();
626  tetgen_in.facetlist = new tetgenio::facet[tetgen_in.numberoffacets];
627  tetgen_in.facetmarkerlist = new int[tetgen_in.numberoffacets];
628  memset(tetgen_in.facetmarkerlist,(int)0,tetgen_in.numberoffacets*sizeof(int));
629  int iHoleInList;      // #Hole for this Facette
630  int iFirstItemPt;
631  // floor
632  iHoleInList = 0;
633  tetgen_in.facetmarkerlist[iFacInList] = 0;
634  FacWall = tetgen_in.facetlist+(iFacInList++);
635  SetFacette(FacWall, 0, 0, 2, 6, 4,_Equip.size()); // Set Facet of wall
636
637  for (iItem=0;iItem<_Equip.size();iItem++) {
638    iFirstItemPt = iPtInList;
639    // Add Points of Item
640    SetPoint(tetgen_in.pointlist,iPtInList++,_Equip[iItem].min_[0],_Equip[iItem].min_[1],_Equip[iItem].min_[2]);
641    SetPoint(tetgen_in.pointlist,iPtInList++,_Equip[iItem].min_[0],_Equip[iItem].min_[1],_Equip[iItem].max_[2]);
642    SetPoint(tetgen_in.pointlist,iPtInList++,_Equip[iItem].min_[0],_Equip[iItem].max_[1],_Equip[iItem].min_[2]);
643    SetPoint(tetgen_in.pointlist,iPtInList++,_Equip[iItem].min_[0],_Equip[iItem].max_[1],_Equip[iItem].max_[2]);
644    SetPoint(tetgen_in.pointlist,iPtInList++,_Equip[iItem].max_[0],_Equip[iItem].min_[1],_Equip[iItem].min_[2]);
645    SetPoint(tetgen_in.pointlist,iPtInList++,_Equip[iItem].max_[0],_Equip[iItem].min_[1],_Equip[iItem].max_[2]);
646    SetPoint(tetgen_in.pointlist,iPtInList++,_Equip[iItem].max_[0],_Equip[iItem].max_[1],_Equip[iItem].min_[2]);
647    SetPoint(tetgen_in.pointlist,iPtInList++,_Equip[iItem].max_[0],_Equip[iItem].max_[1],_Equip[iItem].max_[2]);
648
649    // Set Hole in Wall Facet
650    p = FacWall->polygonlist+1+iHoleInList;
651    p->numberofvertices = 4;
652    p->vertexlist = new int[p->numberofvertices];
653    p->vertexlist[0]=iFirstItemPt+0;
654    p->vertexlist[1]=iFirstItemPt+2;
655    p->vertexlist[2]=iFirstItemPt+6;
656    p->vertexlist[3]=iFirstItemPt+4;
657    FacWall->holelist[3*iHoleInList+0] = 0.5*(_Equip[iItem].min_[0]+_Equip[iItem].max_[0]);
658    FacWall->holelist[3*iHoleInList+1] = 0.5*(_Equip[iItem].min_[1]+_Equip[iItem].max_[1]);
659    FacWall->holelist[3*iHoleInList+2] = 0.0;
660    iHoleInList++;
661
662    // Add inner Facets (Equipment Boundary)
663    SetFacette(tetgen_in.facetlist+(iFacInList++), iFirstItemPt, 1, 3, 7, 5);
664    SetFacette(tetgen_in.facetlist+(iFacInList++), iFirstItemPt, 0, 1, 3, 2);
665    SetFacette(tetgen_in.facetlist+(iFacInList++), iFirstItemPt, 4, 5, 7, 6);
666    SetFacette(tetgen_in.facetlist+(iFacInList++), iFirstItemPt, 0, 1, 5, 4);
667    SetFacette(tetgen_in.facetlist+(iFacInList++), iFirstItemPt, 2, 3, 7, 6);
668  }
669
670  // ceiling
671  iHoleInList = 0;
672  tetgen_in.facetmarkerlist[iFacInList] = 0;
673  FacWall = tetgen_in.facetlist+(iFacInList++);
674  tetgenio::init(FacWall);
675  FacWall->numberofholes=0;
676  FacWall->numberofpolygons=1+FacWall->numberofholes;
677  FacWall->polygonlist = new tetgenio::polygon[FacWall->numberofpolygons];
678  p = FacWall->polygonlist;
679  tetgenio::init(p);
680  p->numberofvertices = 4;
681  p->vertexlist = new int[p->numberofvertices];
682  p->vertexlist[0]=1;
683  p->vertexlist[1]=3;
684  p->vertexlist[2]=7;
685  p->vertexlist[3]=5;
686
687  Item* pItem;
688
689  // x0==0
691  int iWall(0);
692  iHoleInList = 0;
693  tetgen_in.facetmarkerlist[iFacInList] = 0;
694  FacWall = tetgen_in.facetlist+(iFacInList++);
695  SetFacette(FacWall, 0, 0, 1, 3, 2,NumOfWallItems[iWall]);
697  ppItem=_ItemAtWall[iWall].begin();
698  for (;ppItem!=_ItemAtWall[iWall].end();ppItem++) {
699    pItem = *ppItem;
700    iFirstItemPt = iPtInList;
701    SetPoint(tetgen_in.pointlist,iPtInList++,pItem->min_[0],pItem->min_[1],pItem->min_[2]);
702    SetPoint(tetgen_in.pointlist,iPtInList++,pItem->min_[0],pItem->min_[1],pItem->max_[2]);
703    SetPoint(tetgen_in.pointlist,iPtInList++,pItem->min_[0],pItem->max_[1],pItem->min_[2]);
704    SetPoint(tetgen_in.pointlist,iPtInList++,pItem->min_[0],pItem->max_[1],pItem->max_[2]);
705    p = FacWall->polygonlist+1+iHoleInList;
706    p->numberofvertices = 4;
707    p->vertexlist = new int[p->numberofvertices];
708    p->vertexlist[0]=iFirstItemPt+0;
709    p->vertexlist[1]=iFirstItemPt+1;
710    p->vertexlist[2]=iFirstItemPt+3;
711    p->vertexlist[3]=iFirstItemPt+2;
712    FacWall->holelist[3*iHoleInList+0] = 0.0;
713    FacWall->holelist[3*iHoleInList+1] = 0.5*(pItem->min_[1]+pItem->max_[1]);
714    FacWall->holelist[3*iHoleInList+2] = 0.5*(pItem->min_[2]+pItem->max_[2]);
715    iHoleInList++;
716    tetgen_in.facetmarkerlist[iFacInList] = iFacMakerUsed++;
717    SetFacette(tetgen_in.facetlist+(iFacInList++), iFirstItemPt, 0, 1, 3, 2);
718  }
719
720  // x0==_RoomSize[0]
721  iWall++;
722  iHoleInList = 0;
723  tetgen_in.facetmarkerlist[iFacInList] = 0;
724  FacWall = tetgen_in.facetlist+(iFacInList++);
725  SetFacette(FacWall, 0, 4, 5, 7, 6,NumOfWallItems[iWall]);
726  ppItem=_ItemAtWall[iWall].begin();
727  for (;ppItem!=_ItemAtWall[iWall].end();ppItem++) {
728    pItem = *ppItem;
729    iFirstItemPt = iPtInList;
730    SetPoint(tetgen_in.pointlist,iPtInList++,pItem->max_[0],pItem->min_[1],pItem->min_[2]);
731    SetPoint(tetgen_in.pointlist,iPtInList++,pItem->max_[0],pItem->min_[1],pItem->max_[2]);
732    SetPoint(tetgen_in.pointlist,iPtInList++,pItem->max_[0],pItem->max_[1],pItem->min_[2]);
733    SetPoint(tetgen_in.pointlist,iPtInList++,pItem->max_[0],pItem->max_[1],pItem->max_[2]);
734    p = FacWall->polygonlist+1+iHoleInList;
735    p->numberofvertices = 4;
736    p->vertexlist = new int[p->numberofvertices];
737    p->vertexlist[0]=iFirstItemPt+0;
738    p->vertexlist[1]=iFirstItemPt+1;
739    p->vertexlist[2]=iFirstItemPt+3;
740    p->vertexlist[3]=iFirstItemPt+2;
741    FacWall->holelist[3*iHoleInList+0] = _RoomSize[0];
742    FacWall->holelist[3*iHoleInList+1] = 0.5*(pItem->min_[1]+pItem->max_[1]);
743    FacWall->holelist[3*iHoleInList+2] = 0.5*(pItem->min_[2]+pItem->max_[2]);
744    iHoleInList++;
745    tetgen_in.facetmarkerlist[iFacInList] = iFacMakerUsed++;
746    SetFacette(tetgen_in.facetlist+(iFacInList++), iFirstItemPt, 0, 1, 3, 2);
747  }
748
749  // x1==0
750  iWall++;
751  iHoleInList = 0;
752  tetgen_in.facetmarkerlist[iFacInList] = 0;
753  FacWall = tetgen_in.facetlist+(iFacInList++);
754  SetFacette(FacWall, 0, 0, 1, 5, 4, NumOfWallItems[iWall]);
755  ppItem=_ItemAtWall[iWall].begin();
756  for (;ppItem!=_ItemAtWall[iWall].end();ppItem++) {
757    pItem = *ppItem;
758    iFirstItemPt = iPtInList;
759    SetPoint(tetgen_in.pointlist,iPtInList++,pItem->min_[0],pItem->min_[1],pItem->min_[2]);
760    SetPoint(tetgen_in.pointlist,iPtInList++,pItem->min_[0],pItem->min_[1],pItem->max_[2]);
761    SetPoint(tetgen_in.pointlist,iPtInList++,pItem->max_[0],pItem->min_[1],pItem->min_[2]);
762    SetPoint(tetgen_in.pointlist,iPtInList++,pItem->max_[0],pItem->min_[1],pItem->max_[2]);
763    p = FacWall->polygonlist+1+iHoleInList;
764    p->numberofvertices = 4;
765    p->vertexlist = new int[p->numberofvertices];
766    p->vertexlist[0]=iFirstItemPt+0;
767    p->vertexlist[1]=iFirstItemPt+1;
768    p->vertexlist[2]=iFirstItemPt+3;
769    p->vertexlist[3]=iFirstItemPt+2;
770    FacWall->holelist[3*iHoleInList+0] = 0.5*(pItem->min_[0]+pItem->max_[0]);
771    FacWall->holelist[3*iHoleInList+1] = 0.0;
772    FacWall->holelist[3*iHoleInList+2] = 0.5*(pItem->min_[2]+pItem->max_[2]);
773    iHoleInList++;
774    tetgen_in.facetmarkerlist[iFacInList] = iFacMakerUsed++;
775    SetFacette(tetgen_in.facetlist+(iFacInList++), iFirstItemPt, 0, 1, 3, 2);
776  }
777
778  // x1==_RoomSize[1]
779  iWall++;
780  iHoleInList = 0;
781  tetgen_in.facetmarkerlist[iFacInList] = 0;
782  FacWall = tetgen_in.facetlist+(iFacInList++);
783  SetFacette(FacWall, 0, 2, 3, 7, 6, NumOfWallItems[iWall]);
784  ppItem=_ItemAtWall[iWall].begin();
785  for (;ppItem!=_ItemAtWall[iWall].end();ppItem++) {
786    pItem = *ppItem;
787    iFirstItemPt = iPtInList;
788    SetPoint(tetgen_in.pointlist,iPtInList++,pItem->min_[0],pItem->min_[1],pItem->min_[2]);
789    SetPoint(tetgen_in.pointlist,iPtInList++,pItem->min_[0],pItem->min_[1],pItem->max_[2]);
790    SetPoint(tetgen_in.pointlist,iPtInList++,pItem->max_[0],pItem->min_[1],pItem->min_[2]);
791    SetPoint(tetgen_in.pointlist,iPtInList++,pItem->max_[0],pItem->min_[1],pItem->max_[2]);
792    p = FacWall->polygonlist+1+iHoleInList;
793    p->numberofvertices = 4;
794    p->vertexlist = new int[p->numberofvertices];
795    p->vertexlist[0]=iFirstItemPt+0;
796    p->vertexlist[1]=iFirstItemPt+1;
797    p->vertexlist[2]=iFirstItemPt+3;
798    p->vertexlist[3]=iFirstItemPt+2;
799    FacWall->holelist[3*iHoleInList+0] = 0.5*(pItem->min_[0]+pItem->max_[0]);
800    FacWall->holelist[3*iHoleInList+1] = _RoomSize[1];
801    FacWall->holelist[3*iHoleInList+2] = 0.5*(pItem->min_[2]+pItem->max_[2]);
802    iHoleInList++;
803    tetgen_in.facetmarkerlist[iFacInList] = iFacMakerUsed++;
804    SetFacette(tetgen_in.facetlist+(iFacInList++), iFirstItemPt, 0, 1, 3, 2);
805  }
806
807  tetgenio tetgen_out;
808
809  bool bCallExternal=true;
810  if (bCallExternal) {
811    std::ofstream os("Mesh3DRaw.poly",std::ios::out);
812    PrintTetgenMesh(tetgen_in, os);
813    char strBuf[256];
814    sprintf(strBuf,"-nQpqa%f",_h*_h*_h);
815    //sprintf(strBuf,"zpQqa%f",h*h*h);  // tetgen in silent mode
816    char Buf[2014];
817    sprintf(Buf,"tetgen %s Mesh3DRaw.poly",strBuf);
818    std::cout << "calling \"" << Buf << "\"...";
819    std::cout.flush();
820    system(Buf);
821    std::cout << " ... finished" << std::endl;
825  }
826  else {
827    bool PrintMeshingData=true;
828    if (PrintMeshingData) {
829      std::ofstream os("Mesh3DRaw.poly",std::ios::out);
830      PrintTetgenMesh(tetgen_in, os);
831    }
832    tetgenbehavior  tetgen_beh;
833    char strBuf[256];
834    sprintf(strBuf,"-nzQpqa%f",_h*_h*_h);
835    //sprintf(strBuf,"zpQqa%f",h*h*h);  // tetgen in silent mode
836    tetgen_beh.parse_commandline(strBuf);
837    //tetrahedralize(&tetgen_beh, &tetgen_in, &tetgen_out);
838    printf("NOT WORKING\n");
839    exit(-1);
840  }
841
842  Tetgen2Mesh(tetgen_out, p_mesh);
843  p_mesh->prepare_for_use();
844// TODO: Clearify what's going on here, why program termination?
845//  tetgen_in.deinitialize();
846//  tetgen_out.deinitialize();
847}
848
849void ProblemGeometry::CreateMesh2D(libMesh::UnstructuredMesh* p_mesh)
850{
851  assert(GetDim()==2);
852  int dim = GetDim();
853  libMesh::Triangle::triangulateio tri_in, tri_out;
854  libMesh::Triangle::init(tri_in);
855  libMesh::Triangle::init(tri_out);
856
857  tri_in.numberofpoints = 4+4*_Equip.size()+ 2*(_AC.size()+_Exh.size());
858  tri_in.pointlist = static_cast<REAL*>(std::malloc(2*tri_in.numberofpoints*sizeof(REAL)));
859  int iPtInList = 0;
860  tri_in.numberofsegments = 4+4*_Equip.size()+2*(_AC.size()+_Exh.size());
861  tri_in.segmentlist = static_cast<int*> (std::malloc(tri_in.numberofsegments*2*sizeof(int)));
862  int iSegInList = 0;
863
864  ////////////////////////
865  // Outer Boundary (including AC's, Exhausts,
866  // Roundtrip: x1=0 -> x2=1 -> x1=1 -> x2=0
867
868  double epsilon = 1e-8;
869
870  std::vector<double> BdPts;
871  BdPts.reserve(_AC.size()+_Exh.size());
872  int iItem;
873  int iPt;
874  std::vector<CompareItem> BdItems;
875  std::vector<CompareItem>::iterator pSortItem;
876  ValidateItemAtWall();
877
878  BdItems.reserve(_AC.size()+_Exh.size());
879  std::list<Item*>::iterator ppItem;
880
881  int iWall(0);
882  CompareItem item;
883  double start, end;
884
885  tri_in.pointlist[2*iPtInList] = 0.0;
886  tri_in.pointlist[2*iPtInList+1] = 0.0;
887  ++iPtInList;
888
889  // Outer Boundary: x1=0, (0,0)->(0,Max);
890  iWall=0;
891  start=0.0;
892  end=_RoomSize[1];
893  if (_ItemAtWall[iWall].empty()) {
894    tri_in.pointlist[2*iPtInList] = 0.0;
895    tri_in.pointlist[2*iPtInList+1] = end;
896    tri_in.segmentlist[2*iSegInList] = iPtInList-1;
897    tri_in.segmentlist[2*iSegInList+1] = iPtInList;
898    ++iPtInList;
899    ++iSegInList;
900  }
901  else {
902    BdItems.clear();
903    ppItem=_ItemAtWall[iWall].begin();
904    for (;ppItem!=_ItemAtWall[iWall].end();++ppItem) {
905      item.pItem = *ppItem;
906      item.SortValue = (*ppItem)->min_[1];
907      BdItems.push_back(item);
908    }
909    sort(BdItems.begin(), BdItems.end(),CompareItems);
910    for (iItem=0;iItem<BdItems.size();++iItem) {
911      tri_in.pointlist[2*iPtInList] = 0.0;
912      tri_in.pointlist[2*iPtInList+1] = BdItems[iItem].pItem->min_[1];
913      tri_in.segmentlist[2*iSegInList] = iPtInList-1;
914      tri_in.segmentlist[2*iSegInList+1] = iPtInList;
915      ++iPtInList;
916      ++iSegInList;
917      tri_in.pointlist[2*iPtInList] = 0.0;
918      tri_in.pointlist[2*iPtInList+1] = BdItems[iItem].pItem->max_[1];
919      tri_in.segmentlist[2*iSegInList] = iPtInList-1;
920      tri_in.segmentlist[2*iSegInList+1] = iPtInList;
921      ++iPtInList;
922      ++iSegInList;
923    }
924    if (fabs(BdItems.back().pItem->max_[1]-end)>epsilon) {
925      tri_in.pointlist[2*iPtInList] = 0.0;
926      tri_in.pointlist[2*iPtInList+1] = end;
927      tri_in.segmentlist[2*iSegInList] = iPtInList-1;
928      tri_in.segmentlist[2*iSegInList+1] = iPtInList;
929      ++iPtInList;
930      ++iSegInList;
931    }
932  }
933
934  // Outer Boundary: x2=Max, (0,Max)->(Max,Max);
935  iWall = 3;
936  start=0.0;
937  end=_RoomSize[0];
938  if (_ItemAtWall[iWall].empty()) {
939    tri_in.pointlist[2*iPtInList] = end;
940    tri_in.pointlist[2*iPtInList+1] = _RoomSize[1];
941    tri_in.segmentlist[2*iSegInList] = iPtInList-1;
942    tri_in.segmentlist[2*iSegInList+1] = iPtInList;
943    ++iPtInList;
944    ++iSegInList;
945  }
946  else {
947    BdItems.clear();
948    ppItem=_ItemAtWall[iWall].begin();
949    for (;ppItem!=_ItemAtWall[iWall].end();++ppItem) {
950      item.pItem = *ppItem;
951      item.SortValue = (*ppItem)->min_[0];
952      BdItems.push_back(item);
953    }
954    sort(BdItems.begin(), BdItems.end(),CompareItems);
955    for (iItem=0;iItem<BdItems.size();++iItem) {
956      tri_in.pointlist[2*iPtInList] = BdItems[iItem].pItem->min_[0];
957      tri_in.pointlist[2*iPtInList+1] = _RoomSize[1];
958      tri_in.segmentlist[2*iSegInList] = iPtInList-1;
959      tri_in.segmentlist[2*iSegInList+1] = iPtInList;
960      ++iPtInList;
961      ++iSegInList;
962      tri_in.pointlist[2*iPtInList] = BdItems[iItem].pItem->max_[0];
963      tri_in.pointlist[2*iPtInList+1] = _RoomSize[1];
964      tri_in.segmentlist[2*iSegInList] = iPtInList-1;
965      tri_in.segmentlist[2*iSegInList+1] = iPtInList;
966      ++iPtInList;
967      ++iSegInList;
968    }
969    if (fabs(BdItems.back().pItem->max_[0]-end)>epsilon) {
970      tri_in.pointlist[2*iPtInList] = end;
971      tri_in.pointlist[2*iPtInList+1] = _RoomSize[1];
972      tri_in.segmentlist[2*iSegInList] = iPtInList-1;
973      tri_in.segmentlist[2*iSegInList+1] = iPtInList;
974      ++iPtInList;
975      ++iSegInList;
976    }
977  }
978  // Outer Boundary: x1=Max, (Max,Max)->(Max,0);
979  iWall = 1;
980  start=_RoomSize[1];
981  end=0.0;
982  if (_ItemAtWall[iWall].empty()) {
983    tri_in.pointlist[2*iPtInList] = _RoomSize[0];
984    tri_in.pointlist[2*iPtInList+1] = end;
985    tri_in.segmentlist[2*iSegInList] = iPtInList-1;
986    tri_in.segmentlist[2*iSegInList+1] = iPtInList;
987    ++iPtInList;
988    ++iSegInList;
989  }
990  else {
991    BdItems.clear();
992    ppItem=_ItemAtWall[iWall].begin();
993    for (;ppItem!=_ItemAtWall[iWall].end();++ppItem) {
994      item.pItem = *ppItem;
995      item.SortValue = (*ppItem)->min_[1];
996      BdItems.push_back(item);
997    }
998    sort(BdItems.begin(), BdItems.end(),CompareItems);
999    for (iItem=BdItems.size()-1;iItem>=0;--iItem) {
1000      tri_in.pointlist[2*iPtInList] = _RoomSize[0];
1001      tri_in.pointlist[2*iPtInList+1] = BdItems[iItem].pItem->max_[1];
1002      tri_in.segmentlist[2*iSegInList] = iPtInList-1;
1003      tri_in.segmentlist[2*iSegInList+1] = iPtInList;
1004      ++iPtInList;
1005      ++iSegInList;
1006      tri_in.pointlist[2*iPtInList] = _RoomSize[0];
1007      tri_in.pointlist[2*iPtInList+1] = BdItems[iItem].pItem->min_[1];
1008      tri_in.segmentlist[2*iSegInList] = iPtInList-1;
1009      tri_in.segmentlist[2*iSegInList+1] = iPtInList;
1010      ++iPtInList;
1011      ++iSegInList;
1012    }
1013    if (fabs(BdItems.back().pItem->max_[1]-end)>epsilon) {
1014      tri_in.pointlist[2*iPtInList] = _RoomSize[0];
1015      tri_in.pointlist[2*iPtInList+1] = end;
1016      tri_in.segmentlist[2*iSegInList] = iPtInList-1;
1017      tri_in.segmentlist[2*iSegInList+1] = iPtInList;
1018      ++iPtInList;
1019      ++iSegInList;
1020    }
1021  }
1022
1023  // Outer Boundary: x2=0, (Max,0)->(0,0);
1024  iWall = 2;
1025  start=_RoomSize[0];
1026  end=0.0;
1027  if (_ItemAtWall[iWall].empty()) {
1028    tri_in.segmentlist[2*iSegInList] = iPtInList-1;
1029    tri_in.segmentlist[2*iSegInList+1] = 0;
1030    ++iSegInList;
1031  }
1032  else {
1033    BdItems.clear();
1034    ppItem=_ItemAtWall[iWall].begin();
1035    for (;ppItem!=_ItemAtWall[iWall].end();++ppItem) {
1036      item.pItem = *ppItem;
1037      item.SortValue = (*ppItem)->min_[0];
1038      BdItems.push_back(item);
1039    }
1040    sort(BdItems.begin(), BdItems.end(),CompareItems);
1041    for (iItem=BdItems.size()-1;iItem>=0;--iItem) {
1042      tri_in.pointlist[2*iPtInList] = BdItems[iItem].pItem->max_[0];
1043      tri_in.pointlist[2*iPtInList+1] = 0.0;
1044      tri_in.segmentlist[2*iSegInList] = iPtInList-1;
1045      tri_in.segmentlist[2*iSegInList+1] = iPtInList;
1046      ++iPtInList;
1047      ++iSegInList;
1048      tri_in.pointlist[2*iPtInList] = BdItems[iItem].pItem->min_[0];
1049      tri_in.pointlist[2*iPtInList+1] = 0.0;
1050      tri_in.segmentlist[2*iSegInList] = iPtInList-1;
1051      tri_in.segmentlist[2*iSegInList+1] = iPtInList;
1052      ++iPtInList;
1053      ++iSegInList;
1054    }
1055    if (fabs(BdItems.back().pItem->max_[0]-end)>epsilon) {
1056      tri_in.segmentlist[2*iSegInList] = iPtInList-1;
1057      tri_in.segmentlist[2*iSegInList+1] = 0;
1058      ++iSegInList;
1059    }
1060    else {
1061      tri_in.segmentlist[2*iSegInList+1] = 0;
1062    }
1063  }
1064  ///////////////////////////////////////////
1066  tri_in.numberofholes = _Equip.size();
1067  tri_in.holelist = static_cast<REAL*>(std::malloc(2*tri_in.numberofholes*sizeof(REAL)));
1068  int iHoleInList = 0;
1069  Point Center;
1070  for (iItem=0; iItem<_Equip.size();iItem++) {
1071    // Add hole boundary points and segments
1072    tri_in.pointlist[2*iPtInList+0] = _Equip[iItem].min_[0];
1073    tri_in.pointlist[2*iPtInList+1] = _Equip[iItem].min_[1];
1074    tri_in.pointlist[2*iPtInList+2] = _Equip[iItem].min_[0];
1075    tri_in.pointlist[2*iPtInList+3] = _Equip[iItem].max_[1];
1076    tri_in.pointlist[2*iPtInList+4] = _Equip[iItem].max_[0];
1077    tri_in.pointlist[2*iPtInList+5] = _Equip[iItem].min_[1];
1078    tri_in.pointlist[2*iPtInList+6] = _Equip[iItem].max_[0];
1079    tri_in.pointlist[2*iPtInList+7] = _Equip[iItem].max_[1];
1080    tri_in.segmentlist[2*iSegInList+0] = iPtInList+0;
1081    tri_in.segmentlist[2*iSegInList+1] = iPtInList+1;
1082    tri_in.segmentlist[2*iSegInList+2] = iPtInList+1;
1083    tri_in.segmentlist[2*iSegInList+3] = iPtInList+3;
1084    tri_in.segmentlist[2*iSegInList+4] = iPtInList+3;
1085    tri_in.segmentlist[2*iSegInList+5] = iPtInList+2;
1086    tri_in.segmentlist[2*iSegInList+6] = iPtInList+2;
1087    tri_in.segmentlist[2*iSegInList+7] = iPtInList+0;
1088    iPtInList+=4;
1089    iSegInList+=4;
1090    tri_in.holelist[2*iHoleInList+0] = 0.5*(_Equip[iItem].min_[0]+_Equip[iItem].max_[0]);
1091    tri_in.holelist[2*iHoleInList+1] = 0.5*(_Equip[iItem].min_[1]+_Equip[iItem].max_[1]);
1092    iHoleInList++;
1093  }
1094
1095  bool PrintMeshingData=false;
1096  if (PrintMeshingData) {
1097    std::ofstream file("Mesh2DIn.poly");
1098    PrintTriangleMesh(tri_in, file);
1099  }
1100  char strBuf[256];
1101  sprintf(strBuf,"zQnqpa%f",_h*_h);
1102
1103  triangulate(strBuf, &tri_in, &tri_out, NULL);
1104  //libMesh::Triangle::copy_tri_to_mesh(tri_out,*p_mesh,TRI3);
1105  Triangle2Mesh(tri_out,p_mesh);
1106  p_mesh->prepare_for_use();
1107  libMesh::Triangle::destroy(tri_in,Triangle::INPUT);   // deletes also the memory
1108  libMesh::Triangle::destroy(tri_out,Triangle::OUTPUT);
1109}
1110
1112{
1113  char Buf[1024];
1114  double Vals[4];
1115  int n_nodes, dim;
1117
1118  std::ifstream f(str.c_str(),std::ios::in);
1119  f.getline(Buf,1024);
1120
1121  read = sscanf(Buf,"%d %d %lf %lf", &n_nodes,&dim,Vals,Vals+1);
1124    str += Buf;
1125    throw std::runtime_error(str);
1126  }
1127  if (dim!=3) {
1128    std::string str("wrong dimension in node file:");
1129    str += Buf;
1130    throw std::runtime_error(str);
1131  }
1132
1133  tet->numberofpoints = n_nodes;
1134  tet->pointlist = new double[3*tet->numberofpoints];
1135  int i_node;
1136  for (int iPt=0;iPt<tet->numberofpoints;iPt++) {
1137    f.getline(Buf,1024);
1138    read = sscanf(Buf,"%d %lf %lf %lf", &i_node,Vals,Vals+1,Vals+2);
1139    if ( read!=4 ) {
1140      std::string str("Can't read node from line:");
1141      str += Buf;
1142      throw std::runtime_error(str);
1143    }
1144    else {
1145      assert(iPt==i_node-1);
1146      tet->pointlist[3*iPt+0] = Vals[0];
1147      tet->pointlist[3*iPt+1] = Vals[1];
1148      tet->pointlist[3*iPt+2] = Vals[2];
1149    }
1150  }
1151}
1152
1154{
1155  char Buf[1024];
1156  int Vals[5];
1157  int n_elems, n_PtsPerTet;
1159
1160  std::ifstream f(str.c_str(),std::ios::in);
1161  f.getline(Buf,1024);
1162
1163  read = sscanf(Buf,"%d %d %d %d", &n_elems,&n_PtsPerTet,Vals,Vals+1);
1166    str += Buf;
1167    throw std::runtime_error(str);
1168  }
1169  if (n_PtsPerTet!=4) {
1170    std::string str("Wrong number of points per tetrahedron:");
1171    str += Buf;
1172    throw std::runtime_error(str);
1173  }
1174
1175  tet->numberoftetrahedra = n_elems;
1176  tet->tetrahedronlist = new int[4*tet->numberoftetrahedra];
1177  int i_elem;
1178  for (int i_elem=0;i_elem<tet->numberoftetrahedra;i_elem++) {
1179    f.getline(Buf,1024);
1180    read = sscanf(Buf,"%d %d %d %d %d", Vals,Vals+1,Vals+2,Vals+3,Vals+4);
1181    if ( read!=5 ) {
1182      std::string str("Can't read tetrahedron from line:");
1183      str += Buf;
1184      throw std::runtime_error(str);
1185    }
1186    else {
1187      assert(i_elem==Vals[0]-1);
1188      tet->tetrahedronlist[4*i_elem+0] = Vals[1]-1;
1189      tet->tetrahedronlist[4*i_elem+1] = Vals[2]-1;
1190      tet->tetrahedronlist[4*i_elem+2] = Vals[3]-1;
1191      tet->tetrahedronlist[4*i_elem+3] = Vals[4]-1;
1192    }
1193  }
1194}
1195
1197{
1198  char Buf[1024];
1199  int Vals[5];
1200  int n_elems, n_NeighPerTet;
1202
1203  std::ifstream f(str.c_str(),std::ios::in);
1204  f.getline(Buf,1024);
1205
1206  read = sscanf(Buf,"%d %d", &n_elems,&n_NeighPerTet,Vals,Vals+1);
1209    str += Buf;
1210    throw std::runtime_error(str);
1211  }
1212  if (n_NeighPerTet!=4) {
1213    std::string str("Wrong number of neoghnbors per tetrahedron:");
1214    str += Buf;
1215    throw std::runtime_error(str);
1216  }
1217
1218  assert(n_elems==tet->numberoftetrahedra);
1219  tet->neighborlist = new int[4*tet->numberoftetrahedra];
1220  int i_elem;
1221  for (int i_elem=0;i_elem<tet->numberoftetrahedra;i_elem++) {
1222    f.getline(Buf,1024);
1223    read = sscanf(Buf,"%d %d %d %d %d", Vals,Vals+1,Vals+2,Vals+3,Vals+4);
1224    if ( read!=5 ) {
1225      std::string str("Can't read neighbors from line:");
1226      str += Buf;
1227      throw std::runtime_error(str);
1228    }
1229    else {
1230      assert(i_elem==Vals[0]-1);
1231      tet->neighborlist[4*i_elem+0] = Vals[1]<0 ? -1 : (Vals[1]-1);
1232      tet->neighborlist[4*i_elem+1] = Vals[2]<0 ? -1 : (Vals[2]-1);
1233      tet->neighborlist[4*i_elem+2] = Vals[3]<0 ? -1 : (Vals[3]-1);
1234      tet->neighborlist[4*i_elem+3] = Vals[4]<0 ? -1 : (Vals[4]-1);
1235    }
1236  }
1237}
1238
1239// for debugging:
1240void PrintTetgenMesh(const tetgenio& tet, std::ostream& os)
1241{
1242  os.flush();
1243  os << tet.numberofpoints << " " << tet.mesh_dim << std::endl;
1244  os << "# Nodes" << std::endl;
1245  for (int iPt=0;iPt<tet.numberofpoints;iPt++)
1246    os << iPt+1 << " " << tet.pointlist[3*iPt+0] << " " << tet.pointlist[3*iPt+1] << " " << tet.pointlist[3*iPt+2] << std::endl;
1247  os << "# Facets" << std::endl;
1248  tetgenio::facet *f;
1249  tetgenio::polygon *p;
1250
1251  os << tet.numberoffacets << " 1" << std::endl;
1252  for (int iFac=0;iFac<tet.numberoffacets;iFac++) {
1253    f = tet.facetlist+iFac;
1254    os << f->numberofpolygons << " " << f->numberofholes << " " << iFac << std::endl;
1255    for (int iPoly=0; iPoly<f->numberofpolygons; iPoly++) {
1256      p = f->polygonlist+iPoly;
1257      os << p->numberofvertices << " ";
1258      for ( int iPt=0;iPt<p->numberofvertices; iPt++)
1259        os << p->vertexlist[iPt]+1 << " ";
1260      os << std::endl;
1261    }
1262    for (int iPt=0;iPt<f->numberofholes;iPt++) {
1263      os << iPt+1 << " ";
1264      for (int iDim=0;iDim<tet.mesh_dim;iDim++)
1265        os << f->holelist[iPt*tet.mesh_dim+iDim] << " ";
1266      os << std::endl;
1267    }
1268  }
1269
1270  os << "# volume holes" << std::endl << "0" << std::endl << std::endl;
1271  os << "# regions" << std::endl << "0" << std::endl << std::endl;
1272  /*  os << tet.numberoftetrahedra << std::endl;
1273    for(int iEl=0;iEl<tet.numberoftetrahedra;iEl++)
1274      os << iEl << " " << tet.tetrahedronlist[4*iEl+0]+1 << " " << tet.tetrahedronlist[4*iEl+1]+1 << " " << tet.tetrahedronlist[4*iEl+2]+1 << " " << tet.tetrahedronlist[4*iEl+3]+1 << std::endl;
1275    os.flush();
1276  */
1277}
1278
1279//////////////////////////////////////////////////////////////////////////////////////////
1280// 2D
1281void ProblemGeometry::Triangle2Mesh(const libMesh::Triangle::triangulateio& tri, libMesh::UnstructuredMesh* p_mesh)
1282{
1283  const double eps=1e-8;
1284  assert(GetDim()==2);
1285  int dim = GetDim();
1286  p_mesh->set_mesh_dimension(GetDim());
1288  unsigned long int iPt;
1289  p_mesh->reserve_elem(tri.numberoftriangles);
1290  p_mesh->reserve_nodes(tri.numberofpoints);
1291  for (iPt=0;iPt<tri.numberofpoints;iPt++)
1293  // elements
1294  libMesh::Tri3 *pElem;
1295  for (unsigned int iEl=0; iEl<tri.numberoftriangles; iEl++) {
1296    pElem = new libMesh::Tri3;
1297    for (iPt=0;iPt<dim+1;iPt++)
1298      pElem->set_node(iPt) = p_mesh->node_ptr(tri.trianglelist[(dim+1)*iEl+iPt]);
1299    // Finally, add this element to the mesh
1301  }
1302
1303  if (tri.neighborlist) {
1304    for (unsigned int iEl=0; iEl<tri.numberoftriangles; iEl++) {
1305      for (int iNeig=0;iNeig<dim+1;++iNeig) {
1306        if (tri.neighborlist[3*iEl+iNeig]==-1) {
1307          pElem = static_cast<libMesh::Tri3*>(p_mesh->elem(iEl));
1308          int iBdNode1 = tri.trianglelist[(dim+1)*iEl+((iNeig+1)%(dim+1))];
1309          int iBdNode2 = tri.trianglelist[(dim+1)*iEl+((iNeig+2)%(dim+1))];
1310          int sum = iBdNode1+iBdNode2;
1311          for (int iSide=0;iSide<GetDim()+1;iSide++) {
1312            // sum is correct -> sid is correct (holds only, if side has one node less than elem)
1313            if ( pElem->node(pElem->side_nodes_map[iSide][0])+pElem->node(pElem->side_nodes_map[iSide][1]) == sum ) {
1314              std::vector<double> Center;
1315              Center.resize(dim);
1316              for (int iDim=0;iDim<dim;iDim++)
1317                Center[iDim] = 0.5*(tri.pointlist[dim*iBdNode1+iDim] + tri.pointlist[dim*iBdNode2+iDim]);
1318              int BoundaryMarker = GetBoundaryMarker(Center);
1319              assert(BoundaryMarker!=-1);
1321              if ( (fabs(_BoundCond[BoundaryMarker].PhiNeumannCoef)<eps) || (fabs(_BoundCond[BoundaryMarker].TNeumannCoef)<eps) ) {
1324              }
1325            }
1326          }
1327        }
1328      }
1329    }
1330  }
1331  p_mesh->prepare_for_use(false);
1332}
1333
1334void PrintTriangleMesh(const libMesh::Triangle::triangulateio& tri, std::ostream& os)
1335{
1336  os.flush();
1337  os << "# Nodes:" << std::endl;
1338  os << tri.numberofpoints << " 2 0 0" << std::endl;
1339  for (int iPt=0;iPt<tri.numberofpoints; ++iPt)
1340    os << iPt << " " << tri.pointlist[2*iPt] << " " << tri.pointlist[2*iPt+1] << std::endl;
1341  os << std::endl;
1342  os << "# Edges:" << std::endl;
1343  os << tri.numberofsegments << " " << (int)(tri.segmentmarkerlist!=NULL) << std::endl;
1344  for (int iEdge=0;iEdge<tri.numberofsegments;iEdge++) {
1345    os << iEdge << " " << tri.segmentlist[2*iEdge] << " " << tri.segmentlist[2*iEdge+1];
1346    if (tri.segmentmarkerlist!=NULL)
1347      os << " " << tri.segmentmarkerlist[iEdge];
1348    os << std::endl;
1349  }
1350  os << std::endl;
1351  os << "# Holes:" << std::endl;
1352  os << tri.numberofholes << " 0" << std::endl;
1353  for (int iHole=0;iHole<tri.numberofholes;iHole++)
1354    os << iHole << " " << tri.holelist[2*iHole] << " " << tri.holelist[2*iHole+1] << std::endl;
1355  os << std::endl;
1356  os << "# Region attrib., area constr.:" << std::endl;
1357  os << "0";
1358}
1359
1360std::ostream& operator << (std::ostream& os, const BoundaryCondition& BC)
1361{
1362  os << "BC: " << BC.PhiDiricheltCoef << ", " << BC.PhiNeumannCoef << ", " << BC.PhiRhs << ", " << BC.TDiricheltCoef << ", " << BC.TNeumannCoef << ", " << BC.TRhs << std::endl;
1363  return os;
1364}
1365
1366std::ostream& operator << (std::ostream& os, const std::vector<BoundaryCondition>& BCs)
1367{
1368  for (int i=0;i<BCs.size(); i++) os << BCs[i];
1369  return os;
1370}
Note: See TracBrowser for help on using the repository browser.