Source Code for Module meshpy.common

  1 -class _Table: 
  2 -    def __init__(self): 
  3          self.Rows = [] 
  4   
  5 -    def add_row(self, row): 
  6          self.Rows.append([str(i) for i in row]) 
  7   
  8 -    def __str__(self): 
  9          columns = len(self.Rows[0]) 
 10          col_widths = [max(len(row[i]) for row in self.Rows) 
 11                        for i in range(columns)] 
 12   
 13          lines = [ 
 14              " ".join([cell.ljust(col_width) 
 15                        for cell, col_width in zip(row, col_widths)]) 
 16              for row in self.Rows] 
 17          return "\n".join(lines) 
 18   
 19   
 20   
 21   
 22 -def _linebreak_list(list, per_line=10, pad=None): 
 23      def format(s): 
 24          if pad is None: 
 25              return str(s) 
 26          else: 
 27              return str(s).rjust(pad) 
 28   
 29      result = "" 
 30      while len(list) > per_line: 
 31          result += " ".join(format(l) for l in list[:per_line]) + "\n" 
 32          list = list[per_line:] 
 33      return result + " ".join(format(l) for l in list) 
 34           
 35   
 36   
 37 -class MeshInfoBase: 
 38      @property 
 39 -    def face_vertex_indices_to_face_marker(self): 
 40          try: 
 41              return self._fvi2fm 
 42          except AttributeError: 
 43              result = {} 
 44   
 45              for i, face in enumerate(self.faces): 
 46                  result[frozenset(face)] = self.face_markers[i] 
 47   
 48              self._fvi2fm = result 
 49              return result 
 50   
 51   
 52   
 53   
 54   
 55 -    def set_points(self, points, point_markers=None): 
 56          if point_markers is not None: 
 57              assert len(point_markers) == len(point_markers) 
 58   
 59          self.points.resize(len(points)) 
 60   
 61          for i, pt in enumerate(points): 
 62              self.points[i] = pt 
 63         
 64          if point_markers is not None: 
 65              for i, mark in enumerate(point_markers): 
 66                  self.point_markers[i] = mark 
 67   
 68   
 69   
 70   
 71   
 72 -    def set_holes(self, hole_starts): 
 73          self.holes.resize(len(hole_starts)) 
 74          for i, hole in enumerate(hole_starts): 
 75              self.holes[i] = hole 
 76   
 77   
 78   
 79   
 80 -    def write_neu(self, outfile, bc={}, periodicity=None, description="MeshPy Output"): 
 81          """Write the mesh out in (an approximation to) Gambit neutral mesh format. 
 82           
 83          outfile is a file-like object opened for writing. 
 84   
 85          bc is a dictionary mapping single face markers (or frozensets of them) 
 86          to a tuple (bc_name, bc_code). 
 87   
 88          periodicity is either a tuple (face_marker, (px,py,..)) giving the  
 89          face marker of the periodic boundary and the period in each coordinate 
 90          direction (0 if none) or the value None for no periodicity. 
 91          """ 
 92   
 93          from meshpy import version 
 94          from datetime import datetime 
 95   
 96          # header -------------------------------------------------------------- 
 97          outfile.write("CONTROL INFO 2.1.2\n") 
 98          outfile.write("** GAMBIT NEUTRAL FILE\n") 
 99          outfile.write("%s\n" % description) 
100          outfile.write("PROGRAM: MeshPy VERSION: %s\n" % version) 
101          outfile.write("%s\n" % datetime.now().ctime()) 
102           
103          bc_markers = bc.keys() 
104          if periodicity: 
105              periodic_marker, periods = periodicity 
106              bc_markers.append(periodic_marker) 
107               
108          assert len(self.points) 
109   
110          dim = len(self.points[0]) 
111          data = ( 
112                  ("NUMNP", len(self.points)), 
113                  ("NELEM", len(self.elements)), 
114                  ("NGRPS", 1), 
115                  ("NBSETS", len(bc_markers)), 
116                  ("NDFCD", dim), 
117                  ("NDFVL", dim), 
118                  ) 
119   
120          tbl = _Table() 
121          tbl.add_row(key for key, value in data) 
122          tbl.add_row(value for key, value in data) 
123          outfile.write(str(tbl)) 
124          outfile.write("\n") 
125          outfile.write("ENDOFSECTION\n") 
126   
127          # nodes --------------------------------------------------------------- 
128          outfile.write("NODAL COORDINATES 2.1.2\n") 
129          for i, pt in enumerate(self.points): 
130              outfile.write("%d %s\n" %  
131                      (i+1, " ".join(repr(c) for c in pt))) 
132          outfile.write("ENDOFSECTION\n") 
133   
134          # elements ------------------------------------------------------------ 
135          outfile.write("ELEMENTS/CELLS 2.1.2\n") 
136          if dim == 2: 
137              eltype = 3 
138          else: 
139              eltype = 6 
140   
141          for i, el in enumerate(self.elements): 
142              outfile.write("%8d%3d%3d %s\n" %  
143                      (i+1, eltype, len(el),  
144                          "".join("%8d" % (p+1) for p in el))) 
145          outfile.write("ENDOFSECTION\n") 
146   
147          # element groups ------------------------------------------------------ 
148          outfile.write("ELEMENT GROUP 1.3.0\n") 
149          # FIXME 
150          i = 0 
151          grp_elements = range(len(self.elements)) 
152          material = 1 
153          flags = 0 
154          outfile.write("GROUP:%11d ELEMENTS:%11d MATERIAL:%11s NFLAGS: %11d\n" 
155                  % (1, len(grp_elements), repr(material), flags)) 
156          outfile.write(("epsilon: %s\n" % material).rjust(32)) # FIXME 
157          outfile.write("0\n") 
158          outfile.write(_linebreak_list([str(i+1) for i in grp_elements], 
159              pad=8) 
160                  + "\n") 
161          outfile.write("ENDOFSECTION\n") 
162   
163          # boundary conditions ------------------------------------------------- 
164          # build mapping face -> (tet, neu_face_index) 
165          face2el = {} 
166   
167          if dim == 2: 
168              for ti, el in enumerate(self.elements): 
169                  # Sledge++ Users' Guide, figure 4 
170                  faces = [ 
171                          frozenset([el[0], el[1]]), 
172                          frozenset([el[1], el[2]]), 
173                          frozenset([el[2], el[0]]), 
174                          ] 
175                  for fi, face in enumerate(faces): 
176                      face2el.setdefault(face, []).append((ti, fi+1)) 
177   
178          elif dim == 3: 
179              face2el = {} 
180              for ti, el in enumerate(self.elements): 
181                  # Sledge++ Users' Guide, figure 5 
182                  faces = [ 
183                          frozenset([el[1], el[0], el[2]]), 
184                          frozenset([el[0], el[1], el[3]]), 
185                          frozenset([el[1], el[2], el[3]]), 
186                          frozenset([el[2], el[0], el[3]]), 
187                          ] 
188                  for fi, face in enumerate(faces): 
189                      face2el.setdefault(face, []).append((ti, fi+1)) 
190   
191          else: 
192              raise ValueError, "invalid number of dimensions (%d)" % dim 
193   
194          # actually output bc sections 
195          if not self.faces.allocated: 
196              from warnings import warn 
197              warn("no exterior faces in mesh data structure, not writing boundary conditions") 
198          else: 
199              # requires -f option in tetgen, -e in triangle 
200   
201              for bc_marker in bc_markers: 
202                  if isinstance(bc_marker, frozenset): 
203                      face_indices = [i 
204                              for i, face in enumerate(self.faces) 
205                              if self.face_markers[i] in bc_marker] 
206                  else: 
207                      face_indices = [i 
208                              for i, face in enumerate(self.faces) 
209                              if bc_marker == self.face_markers[i]] 
210   
211                  if not face_indices: 
212                      continue 
213   
214                  outfile.write("BOUNDARY CONDITIONS 2.1.2\n") 
215                  if bc_marker in bc: 
216                      # regular BC 
217   
218                      bc_name, bc_code = bc[bc_marker] 
219                      outfile.write("%32s%8d%8d%8d%8d\n"  
220                              % (bc_name,  
221                                  1, # face BC 
222                                  len(face_indices), 
223                                  0, # zero additional values per face, 
224                                  bc_code, 
225                                  ) 
226                              ) 
227                  else: 
228                      # periodic BC 
229   
230                      outfile.write("%s%s%8d%8d%8d\n" 
231                              % ("periodic", " ".join(repr(p) for p in periods), 
232                                  len(face_indices), 
233                                  0, # zero additional values per face, 
234                                  0, 
235                                  ) 
236                              ) 
237   
238                  for i, fi in enumerate(face_indices): 
239                      face_nodes = frozenset(self.faces[fi]) 
240                      adj_el = face2el[face_nodes] 
241                      assert len(adj_el) == 1 
242   
243                      el_index, el_face_number = adj_el[0] 
244   
245                      outfile.write("%10d%5d%5d\n" %  
246                              (el_index+1, eltype, el_face_number)) 
247   
248                  outfile.write("ENDOFSECTION\n") 
249   
250              outfile.close() 
251              # FIXME curved boundaries? 
252              # FIXME proper element group support 
253   
254   
255   
256   
257   
258 -def dump_array(name, array): 
259      print "array %s: %d elements, %d values per element" % (name, len(array), array.unit) 
260   
261      if len(array) == 0 or array.unit == 0: 
262          return 
263   
264      try: 
265          array[0] 
266      except RuntimeError: 
267          print "  not allocated" 
268          return 
269   
270      for i, entry in enumerate(array): 
271          if isinstance(entry, list): 
272              print "  %d: %s" % (i, ",".join(str(sub) for sub in entry)) 
273          else: 
274              print "  %d: %s" % (i, entry) 
275