MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
mesquite.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "mesquite.hpp"
13 
14 #ifdef MFEM_USE_MESQUITE
15 
16 #include "../fem/fem.hpp"
17 #include <iostream>
18 
19 namespace mfem
20 {
21 
22 using namespace std;
23 
25 {
26  if (elementData)
27  {
28  free(elementData);
29  }
30  if (vertexData)
31  {
32  free(vertexData);
33  }
34  if (defaultValue)
35  {
36  free(defaultValue);
37  }
38 }
39 
40 
41 void MesquiteMesh::MeshTags::clear()
42 {
43  for (std::vector<TagData*>::iterator iter = tagList.begin();
44  iter != tagList.end(); ++iter)
45  if (*iter)
46  {
47  delete *iter;
48  }
49 
50  tagList.clear();
51 }
52 
53 size_t MesquiteMesh::MeshTags::size_from_tag_type( Mesh::TagType type )
54 {
55  switch ( type )
56  {
57  case Mesh::BYTE: return 1;
58  case Mesh::BOOL: return sizeof(bool);
59  case Mesh::DOUBLE: return sizeof(double);
60  case Mesh::INT: return sizeof(int);
61  case Mesh::HANDLE: return sizeof(void*);
62  case Mesh::LONG_LONG: return sizeof(long long);
63  default: assert(0); return 0;
64  }
65 }
66 
67 size_t MesquiteMesh::MeshTags::create( const std::string& name,
68  Mesh::TagType type,
69  unsigned length,
70  const void* defval,
71  MsqError& err )
72 {
73  size_t h = handle( name, err );
74  if (h)
75  {
76  MSQ_SETERR(err)(name, MsqError::TAG_ALREADY_EXISTS);
77  return 0;
78  }
79 
80  if (length == 0 || size_from_tag_type(type) == 0)
81  {
82  MSQ_SETERR(err)(MsqError::INVALID_ARG);
83  return 0;
84  }
85 
86  TagData* tag = new TagData( name, type, length );
87  h = tagList.size();
88  tagList.push_back(tag);
89 
90  if (defval)
91  {
92  tag->defaultValue = malloc( tag->desc.size );
93  memcpy( tag->defaultValue, defval, tag->desc.size );
94  }
95 
96  return h+1;
97 }
98 
99 size_t MesquiteMesh::MeshTags::create( const MfemTagDescription& desc,
100  const void* defval,
101  MsqError& err )
102 {
103  size_t h = handle( desc.name.c_str(), err );
104  if (h)
105  {
106  MSQ_SETERR(err)(desc.name.c_str(), MsqError::TAG_ALREADY_EXISTS);
107  return 0;
108  }
109 
110  err.clear();
111  if (desc.size == 0 || (desc.size % size_from_tag_type(desc.type)) != 0)
112  {
113  MSQ_SETERR(err)(MsqError::INVALID_ARG);
114  return 0;
115  }
116 
117  TagData* tag = new TagData( desc );
118  h = tagList.size();
119  tagList.push_back(tag);
120 
121  if (defval)
122  {
123  tag->defaultValue = malloc( tag->desc.size );
124  memcpy( tag->defaultValue, defval, tag->desc.size );
125  }
126 
127  return h+1;
128 }
129 
130 void MesquiteMesh::MeshTags::destroy( size_t tag_index, MsqError& err )
131 {
132  --tag_index;
133  if (tag_index >= tagList.size() || 0 == tagList[tag_index])
134  {
135  MSQ_SETERR(err)(MsqError::TAG_NOT_FOUND);
136  return ;
137  }
138 
139  delete tagList[tag_index];
140  tagList[tag_index] = 0;
141 }
142 
143 size_t MesquiteMesh::MeshTags::handle( const std::string& name,
144  MsqError& err ) const
145 {
146  for (size_t i = 0; i < tagList.size(); ++i)
147  if (tagList[i] && tagList[i]->desc.name == name)
148  {
149  return i+1;
150  }
151 
152  return 0;
153 }
154 
155 const MesquiteMesh::MfemTagDescription& MesquiteMesh::MeshTags::properties(
156  size_t tag_index, MsqError& err ) const
157 {
158  static MfemTagDescription dummy_desc;
159  --tag_index;
160 
161  if (tag_index >= tagList.size() || !tagList[tag_index])
162  {
163  MSQ_SETERR(err)("Invalid tag handle", MsqError::INVALID_ARG);
164  return dummy_desc;
165  }
166 
167  return tagList[tag_index]->desc;
168 }
169 
170 
171 void MesquiteMesh::MeshTags::set_element_data( size_t tag_index,
172  size_t num_indices,
173  const size_t* index_array,
174  const void* values,
175  MsqError& err )
176 {
177  size_t i;
178  char* data;
179  --tag_index;
180  if (tag_index >= tagList.size() || !tagList[tag_index])
181  {
182  MSQ_SETERR(err)("Invalid tag handle", MsqError::INVALID_ARG);
183  return;
184  }
185 
186  TagData* tag = tagList[tag_index];
187 
188  // Get highest element index
189  size_t total = tag->elementCount;
190  for (i = 0; i < num_indices; ++i)
191  if (index_array[i] >= total)
192  {
193  total = index_array[i] + 1;
194  }
195 
196  // If need more space
197  if (total > tag->elementCount)
198  {
199  // allocate more space
200  tag->elementData = realloc( tag->elementData, tag->desc.size * total );
201  // if a default value, initialize new space with it
202  if (tag->defaultValue)
203  {
204  data = ((char*)tag->elementData) + tag->elementCount * tag->desc.size;
205  for (i = tag->elementCount; i < total; ++i)
206  {
207  memcpy( data, tag->defaultValue, tag->desc.size );
208  data += tag->desc.size;
209  }
210  }
211  else
212  {
213  memset( (char*)tag->elementData + tag->elementCount * tag->desc.size, 0,
214  (total - tag->elementCount) * tag->desc.size );
215  }
216  tag->elementCount = total;
217  }
218 
219  // Store passed tag values
220  data = (char*)tag->elementData;
221  const char* iter = (const char*)values;
222  for (i = 0; i < num_indices; ++i)
223  {
224  memcpy( data + index_array[i]*tag->desc.size, iter, tag->desc.size );
225  iter += tag->desc.size;
226  }
227 }
228 
229 void MesquiteMesh::MeshTags::get_element_data( size_t tag_index,
230  size_t num_indices,
231  const size_t* index_array,
232  void* values,
233  MsqError& err ) const
234 {
235  --tag_index;
236  if (tag_index >= tagList.size() || !tagList[tag_index])
237  {
238  MSQ_SETERR(err)("Invalid tag handle", MsqError::INVALID_ARG);
239  return;
240  }
241 
242  TagData* tag = tagList[tag_index];
243 
244  char* iter = (char*)values;
245  const char* data = (const char*)tag->elementData;
246 
247  for (size_t i = 0; i < num_indices; ++i)
248  {
249  const void* ptr;
250  size_t index = index_array[i];
251  if (index >= tag->elementCount)
252  {
253  ptr = tag->defaultValue;
254  if (!ptr)
255  {
256  MSQ_SETERR(err)(MsqError::TAG_NOT_FOUND);
257  return;
258  }
259  }
260  else
261  {
262  ptr = data + index * tag->desc.size;
263  }
264 
265  memcpy( iter, ptr, tag->desc.size );
266  iter += tag->desc.size;
267  }
268 }
269 
270 void MesquiteMesh::MeshTags::set_vertex_data( size_t tag_index,
271  size_t num_indices,
272  const size_t* index_array,
273  const void* values,
274  MsqError& err )
275 {
276  size_t i;
277  char* data;
278  --tag_index;
279  if (tag_index >= tagList.size() || !tagList[tag_index])
280  {
281  MSQ_SETERR(err)("Invalid tag handle", MsqError::INVALID_ARG);
282  return;
283  }
284 
285  TagData* tag = tagList[tag_index];
286 
287  // Get highest element index
288  size_t total = tag->vertexCount;
289  for (i = 0; i < num_indices; ++i)
290  if (index_array[i] >= total)
291  {
292  total = index_array[i] + 1;
293  }
294 
295  // If need more space
296  if (total > tag->vertexCount)
297  {
298  // allocate more space
299  tag->vertexData = realloc( tag->vertexData, tag->desc.size * total );
300  // if a default value, initialize new space with it
301  if (tag->defaultValue)
302  {
303  data = ((char*)tag->vertexData) + tag->vertexCount * tag->desc.size;
304  for (i = tag->vertexCount; i < total; ++i)
305  {
306  memcpy( data, tag->defaultValue, tag->desc.size );
307  data += tag->desc.size;
308  }
309  }
310  else
311  {
312  memset( (char*)tag->vertexData + tag->vertexCount * tag->desc.size, 0,
313  (total - tag->vertexCount) * tag->desc.size );
314  }
315  tag->vertexCount = total;
316  }
317 
318  // Store passed tag values
319  data = (char*)tag->vertexData;
320  const char* iter = (const char*)values;
321  for (i = 0; i < num_indices; ++i)
322  {
323  memcpy( data + index_array[i]*tag->desc.size, iter, tag->desc.size );
324  iter += tag->desc.size;
325  }
326 }
327 
328 void MesquiteMesh::MeshTags::get_vertex_data( size_t tag_index,
329  size_t num_indices,
330  const size_t* index_array,
331  void* values,
332  MsqError& err ) const
333 {
334  --tag_index;
335  if (tag_index >= tagList.size() || !tagList[tag_index])
336  {
337  MSQ_SETERR(err)("Invalid tag handle", MsqError::INVALID_ARG);
338  return;
339  }
340 
341  TagData* tag = tagList[tag_index];
342 
343  char* iter = (char*)values;
344  const char* data = (const char*)tag->vertexData;
345 
346  for (size_t i = 0; i < num_indices; ++i)
347  {
348  const void* ptr;
349  size_t index = index_array[i];
350  if (index >= tag->vertexCount)
351  {
352  ptr = tag->defaultValue;
353  if (!ptr)
354  {
355  MSQ_SETERR(err)(MsqError::TAG_NOT_FOUND);
356  return;
357  }
358  }
359  else
360  {
361  ptr = data + index * tag->desc.size;
362  }
363 
364  memcpy( iter, ptr, tag->desc.size );
365  iter += tag->desc.size;
366  }
367 }
368 
369 bool MesquiteMesh::MeshTags::tag_has_vertex_data( size_t tag_index,
370  MsqError& err )
371 {
372  --tag_index;
373  if (tag_index >= tagList.size() || !tagList[tag_index])
374  {
375  MSQ_SETERR(err)("Invalid tag handle", MsqError::INVALID_ARG);
376  return false;
377  }
378 
379  TagData* tag = tagList[tag_index];
380  return 0 != tag->vertexData || tag->defaultValue;
381 }
382 
383 bool MesquiteMesh::MeshTags::tag_has_element_data( size_t tag_index,
384  MsqError& err )
385 {
386  --tag_index;
387  if (tag_index >= tagList.size() || !tagList[tag_index])
388  {
389  MSQ_SETERR(err)("Invalid tag handle", MsqError::INVALID_ARG);
390  return false;
391  }
392 
393  TagData* tag = tagList[tag_index];
394  return 0 != tag->elementData || tag->defaultValue;
395 }
396 
397 MesquiteMesh::MeshTags::TagIterator MesquiteMesh::MeshTags::tag_begin()
398 {
399  size_t index = 0;
400  while (index < tagList.size() && tagList[index] == NULL)
401  {
402  ++index;
403  }
404  return TagIterator( this, index );
405 }
406 
407 MesquiteMesh::MeshTags::TagIterator
409 {
410  ++index;
411  while (index < tags->tagList.size() && NULL == tags->tagList[index])
412  {
413  ++index;
414  }
415  return TagIterator( tags, index );
416 }
417 
420 {
421  --index;
422  while (index < tags->tagList.size() && NULL == tags->tagList[index])
423  {
424  --index;
425  }
426  return TagIterator( tags, index );
427 }
428 
431 {
432  size_t old = index;
433  ++index;
434  while (index < tags->tagList.size() && NULL == tags->tagList[index])
435  {
436  ++index;
437  }
438  return TagIterator( tags, old );
439 }
440 
443 {
444  size_t old = index;
445  --index;
446  while (index < tags->tagList.size() && NULL == tags->tagList[index])
447  {
448  --index;
449  }
450  return TagIterator( tags, old );
451 }
452 
453 //
454 // MesquiteMesh implementation follows
455 //
457  : myTags( new MeshTags )
458 {
459  mesh = mfem_mesh;
460  nelems = mesh->GetNE();
461  nodes = mesh->GetNodes();
462 
463  if (nodes)
464  {
465  fes = nodes->FESpace();
466  ndofs = fes->GetNDofs();
467  fes->BuildElementToDofTable();
468  const Table *elem_dof = &fes->GetElementToDofTable();
469  dof_elem = new Table;
470  Transpose(*elem_dof, *dof_elem, ndofs);
471  }
472  else
473  {
474  ndofs = mesh->GetNV();
475  dof_elem = mesh->GetVertexToElementTable();
476  }
477 
478  mByte = vector<char>(ndofs);
479  mFixed = vector<bool>(ndofs, false);
480 
481  // By default, flag all boundary nodes as fixed
482  Array<int> bdofs;
483  for (int i = 0; i < mesh->GetNBE(); i++)
484  {
485  if (nodes)
486  {
487  fes->GetBdrElementDofs(i, bdofs);
488  }
489  else
490  {
491  mesh->GetBdrElementVertices(i, bdofs);
492  }
493  for (int j = 0; j < bdofs.Size(); j++)
494  {
495  mFixed[bdofs[j]] = true;
496  }
497  }
498 
499 }
500 
501 int MesquiteMesh::get_geometric_dimension(Mesquite::MsqError &err)
502 {
503  return mesh->Dimension();
504 }
505 
506 void MesquiteMesh::get_all_elements(std::vector<ElementHandle>& elements,
507  Mesquite::MsqError& err)
508 {
509  elements.resize(nelems);
510  for (int i = 0; i < nelems; i++)
511  {
512  elements[i] = (ElementHandle) i;
513  }
514 }
515 
516 void MesquiteMesh::get_all_vertices(std::vector<VertexHandle>& vertices,
517  Mesquite::MsqError& err)
518 {
519  vertices.resize(ndofs);
520  for (int i = 0; i < ndofs; i++)
521  {
522  vertices[i] = (VertexHandle) i;
523  }
524 }
525 
526 void MesquiteMesh::vertices_get_coordinates(const VertexHandle vert_array[],
527  MsqVertex* coordinates,
528  size_t num_vtx,
529  MsqError &err)
530 {
531 
532  const size_t *indices = (const size_t*) vert_array;
533  double coords[3];
534  for (int i = 0; i < num_vtx; i++)
535  {
536  mesh->GetNode(indices[i], coords);
537  coordinates[i].x(coords[0]);
538  coordinates[i].y(coords[1]);
539  if (mesh->Dimension() == 3 )
540  {
541  coordinates[i].z(coords[2]);
542  }
543  else
544  {
545  coordinates[i].z(0.0);
546  }
547  }
548 }
549 
550 void MesquiteMesh::vertex_set_coordinates(VertexHandle vertex,
551  const Vector3D &coordinates,
552  MsqError &err)
553 {
554  double coords[3];
555  coords[0] = coordinates.x();
556  coords[1] = coordinates.y();
557  coords[2] = coordinates.z();
558  mesh->SetNode((size_t) vertex, coords);
559 }
560 
561 void MesquiteMesh::vertex_set_byte(VertexHandle vertex,
562  unsigned char byte,
563  MsqError &err)
564 {
565  size_t index = (size_t) vertex;
566  mByte[index] = byte;
567 }
568 
569 void MesquiteMesh::vertices_set_byte(const VertexHandle *vert_array,
570  const unsigned char *byte_array,
571  size_t array_size,
572  MsqError &err)
573 {
574  const size_t* indices = (const size_t*) vert_array;
575  for (int i = 0; i < array_size; i++)
576  {
577  mByte[indices[i]] = byte_array[i];
578  }
579 }
580 
581 void MesquiteMesh::vertex_get_byte(const VertexHandle vertex,
582  unsigned char *byte,
583  MsqError &err )
584 {
585  *byte = mByte[(const size_t) vertex];
586 }
587 
588 void MesquiteMesh::vertices_get_byte(const VertexHandle *vertex,
589  unsigned char *byte_array,
590  size_t array_size,
591  MsqError &err )
592 {
593  const size_t* indices = (const size_t*) vertex;
594  for (int i = 0; i < array_size; i++)
595  {
596  byte_array[i] = mByte[indices[i]];
597  }
598 }
599 
600 void MesquiteMesh::vertices_get_fixed_flag(const VertexHandle vert_array[],
601  std::vector<bool>& fixed_flag_array,
602  size_t num_vtx,
603  MsqError &err )
604 {
605  fixed_flag_array.resize(num_vtx + 1);
606  const size_t* indices = (const size_t*) vert_array;
607  for (int i = 0; i < num_vtx; i++)
608  {
609  fixed_flag_array[i] = mFixed[indices[i]];
610  }
611 }
612 
613 void MesquiteMesh::vertices_set_fixed_flag(const VertexHandle vert_array[],
614  const std::vector< bool > &fixed_flag_array,
615  size_t num_vtx,
616  MsqError &err )
617 {
618  const size_t* indices = (const size_t*) vert_array;
619  for (int i = 0; i < num_vtx; i++)
620  {
621  mFixed[indices[i]] = fixed_flag_array[i];
622  }
623 }
624 
626  *elem_handles,
627  size_t num_elems,
628  std::vector<VertexHandle>& vert_handles,
629  std::vector<size_t>& offsets,
630  MsqError &err)
631 {
632  const size_t* indices = (const size_t*) elem_handles;
633  vert_handles.clear();
634  offsets.resize(num_elems + 1);
635 
636  Array<int> elem_dofs;
637 
638  for (int i = 0; i < num_elems; i++)
639  {
640  offsets[i] = vert_handles.size();
641  elem = mesh->GetElement(indices[i]);
642  if (nodes)
643  {
644  fes->GetElementDofs(indices[i],elem_dofs);
645  }
646  else
647  {
648  elem->GetVertices(elem_dofs);
649  }
650 
651  for (int j = 0; j < elem_dofs.Size(); j++)
652  {
653  // Ordering of this matters!!!
654  // We are good for triangles, quads and hexes. What about tets?
655  vert_handles.push_back( (VertexHandle) elem_dofs[j] );
656  }
657  }
658  offsets[num_elems] = vert_handles.size();
659 }
660 
662  vertex_array,
663  size_t num_vertex,
664  std::vector<ElementHandle>& elements,
665  std::vector<size_t>& offsets,
666  MsqError& err)
667 {
668  const size_t* indices = (const size_t*) vertex_array;
669  elements.clear();
670  offsets.resize(num_vertex + 1);
671 
672  for (int i = 0; i < num_vertex; i++)
673  {
674  offsets[i] = elements.size();
675  int* vertex_elems = dof_elem->GetRow(indices[i]);
676  for (int j = 0; j < dof_elem->RowSize(indices[i]); j++)
677  {
678  elements.push_back( (ElementHandle) vertex_elems[j] );
679  }
680  }
681  offsets[num_vertex] = elements.size();
682 }
683 
685  *element_handle_array,
686  EntityTopology *element_topologies,
687  size_t num_elements,
688  MsqError &err)
689 {
690  // In MESQUITE:
691  // TRIANGLE = 8
692  // QUADRILATERAL = 9
693  // TETRAHEDRON = 11
694  // HEXAHEDRON = 12
695 
696  // In MFEM:
697  // POINT = 0
698  // SEGMENT = 1
699  // TRIANGLE = 2
700  // QUADRILATERAL = 3
701  // TETRAHEDRON = 4
702  // HEXAHEDRON = 5
703  // BISECTED = 6
704  // QUADRISECTED = 7
705  // OCTASECTED = 8
706 
707  int mfem_to_mesquite[9] = {0,0,8,9,11,12,0,0,0};
708  const size_t* indices = (const size_t*) element_handle_array;
709  for (int i = 0; i < num_elements; i++)
710  {
711  element_topologies[i] = (EntityTopology) mfem_to_mesquite[mesh->GetElementType(
712  indices[i])];
713  }
714 }
715 
717 {
718  MsqPrintError err(mfem::err);
719  delete myTags;
720 
721  delete dof_elem;
722 }
723 
724 void
726 {
727  MsqError err;
728 
729  // create a tag for a single integer value
730  TagHandle attributeTagHandle = tag_create( "material", Mesh::INT, 1, 0, err );
731 
732  int *materialValues = new int[nelems];
733  for (int i=0; i<nelems; i++)
734  {
735  materialValues[i] = mesh->GetAttribute(i);
736  }
737 
738  //
739  // now put these values into an element tag
740  //
741  std::vector<ElementHandle> elements;
742  get_all_elements(elements, err );
743 
744  tag_set_element_data( attributeTagHandle, nelems, arrptr(elements),
745  (const void*)(materialValues), err );
746 
747  delete[] materialValues;
748 }
749 
750 TagHandle MesquiteMesh::tag_create( const std::string& name,
751  TagType type,
752  unsigned length,
753  const void* defval,
754  MsqError& err )
755 {
756 
757  size_t size = MeshTags::size_from_tag_type( type );
758  MfemTagDescription desc( name, type, length*size );
759  size_t index = myTags->create( desc, defval, err ); MSQ_ERRZERO(err);
760  return (TagHandle)index;
761 }
762 
763 void MesquiteMesh::tag_delete( TagHandle handle, MsqError& err )
764 {
765  myTags->destroy( (size_t)handle, err ); MSQ_CHKERR(err);
766 }
767 
768 TagHandle MesquiteMesh::tag_get( const std::string& name, MsqError& err )
769 {
770  size_t index = myTags->handle( name, err ); MSQ_ERRZERO(err);
771  if (!index)
772  {
773  MSQ_SETERR(err)( MsqError::TAG_NOT_FOUND, "could not find tag \"%s\"",
774  name.c_str() );
775  }
776  return (TagHandle)index;
777 }
778 
779 void MesquiteMesh::tag_properties( TagHandle handle,
780  std::string& name,
781  TagType& type,
782  unsigned& length,
783  MsqError& err )
784 {
785  const MfemTagDescription& desc
786  = myTags->properties( (size_t)handle, err ); MSQ_ERRRTN(err);
787 
788  name = desc.name;
789  type = desc.type;
790  length = (unsigned)(desc.size / MeshTags::size_from_tag_type( desc.type ));
791 }
792 
793 
794 void MesquiteMesh::tag_set_element_data( TagHandle handle,
795  size_t num_elems,
796  const ElementHandle* elem_array,
797  const void* values,
798  MsqError& err )
799 {
800  myTags->set_element_data( (size_t)handle,
801  num_elems,
802  (const size_t*)elem_array,
803  values,
804  err ); MSQ_CHKERR(err);
805 }
806 
807 void MesquiteMesh::tag_get_element_data( TagHandle handle,
808  size_t num_elems,
809  const ElementHandle* elem_array,
810  void* values,
811  MsqError& err )
812 {
813  myTags->get_element_data( (size_t)handle,
814  num_elems,
815  (const size_t*)elem_array,
816  values,
817  err ); MSQ_CHKERR(err);
818 }
819 
820 void MesquiteMesh::tag_set_vertex_data( TagHandle handle,
821  size_t num_elems,
822  const VertexHandle* elem_array,
823  const void* values,
824  MsqError& err )
825 {
826  myTags->set_vertex_data( (size_t)handle,
827  num_elems,
828  (const size_t*)elem_array,
829  values,
830  err ); MSQ_CHKERR(err);
831 }
832 
833 void MesquiteMesh::tag_get_vertex_data( TagHandle handle,
834  size_t num_elems,
835  const VertexHandle* elem_array,
836  void* values,
837  MsqError& err )
838 {
839  myTags->get_vertex_data( (size_t)handle,
840  num_elems,
841  (const size_t*)elem_array,
842  values,
843  err ); MSQ_CHKERR(err);
844 }
845 
846 static void BoundaryPreservingOptimization(mfem::MesquiteMesh &mesh)
847 {
848 
849  MsqDebug::enable(1);
850  MsqPrintError err(mfem::err);
851 
852  int pOrder = 2;
853  int mNumInterfaceSmoothIters = 5;
854  bool mFixBndryInInterfaceSmooth = true; //this fixes exterior surface nodes
855  bool project_gradient= false;
856  double cos_crease_angle=0.2;
857 
858  // get all vertices
859  std::vector<Mesquite::Mesh::VertexHandle> vertices;
860  mesh.get_all_vertices(vertices, err);
861  if (MSQ_CHKERR(err)) {mfem::out << err << std::endl; exit(EXIT_FAILURE);}
862  int num_vertices = vertices.size();
863 
864  // get application fixed vertices
865  std::vector<bool> app_fixed(num_vertices);
866 
867  mesh.vertices_get_fixed_flag(&(vertices[0]), app_fixed, num_vertices, err);
868  if (MSQ_CHKERR(err)) {mfem::out << err << std::endl; exit(EXIT_FAILURE);}
869  int num_app_fixed = 0;
870 
871  for (int i = 0; i < num_vertices; i++)
872  {
873  if (app_fixed[i]) { num_app_fixed++; }
874  }
875 
876  mfem::out << "mesh has " << num_vertices << " vertices and " << num_app_fixed <<
877  " are app fixed. ";
878 
879  // create planar domain for interior and assessor queues
880  Mesquite::PlanarDomain geom( PlanarDomain::XY );
881 
882  // tag the underlying mesh material attributes.
883  mesh.tag_attributes();
884 
885  // create boundary domain and find mesh boundary
886  Mesquite::MeshBoundaryDomain2D* mesh_domain = new
887  Mesquite::MeshBoundaryDomain2D( MeshBoundaryDomain2D::XY, 0.0, project_gradient,
888  Mesquite::MeshBoundaryDomain2D::QUADRATIC);
889  mesh_domain->skin_area_mesh(&mesh,cos_crease_angle,"material");
890 
891  std::vector<Mesquite::Mesh::VertexHandle> theBoundaryVertices;
892  mesh_domain->get_boundary_vertices( theBoundaryVertices );
893 
894  int num_boundary_vertices = theBoundaryVertices.size();
895 
896  std::vector<Mesquite::Mesh::VertexHandle> theBoundaryEdges;
897  mesh_domain->get_boundary_edges( theBoundaryEdges );
898 
899  // int num_boundary_edges = theBoundaryEdges.size();
900  std::vector<bool> fixed_flags_boundary(num_boundary_vertices);
901 
902  // get application fixed boundary vertices
903  std::vector<bool> app_fixed_boundary(num_boundary_vertices);
904  mesh.vertices_get_fixed_flag(&(theBoundaryVertices[0]),app_fixed_boundary,
905  num_boundary_vertices, err);
906 
907  int num_app_fixed_boundary = 0;
908  for (int i = 0; i < num_boundary_vertices; i++)
909  {
910  if (app_fixed_boundary[i]) { num_app_fixed_boundary++; }
911  }
912 
913  mfem::out << "mesh has " << num_boundary_vertices << " boundary vertices and "
914  << num_app_fixed_boundary << " are app fixed" << std::endl;
915 
916 
917  // only fix boundary vertices along corners
918  int num_fixed_boundary_flags = 0;
919 
920  std::vector<Mesquite::Mesh::VertexHandle> theCornerVertices;
921  mesh_domain->get_corner_vertices( theCornerVertices );
922 
923 
924  // fix only vertices that are classified as corners
925  for (int i = 0; i < num_boundary_vertices; i++)
926  {
927  if (!mFixBndryInInterfaceSmooth)
928  {
929  fixed_flags_boundary[i] = false;
930  }
931  else
932  {
933  fixed_flags_boundary[i] = app_fixed_boundary[i];
934  }
935 
936  for (int j = 0; j < theCornerVertices.size(); j++)
937  {
938  // printf("theCornerVertices[%d]=%lu\n",j, (size_t)(theCornerVertices[j]));
939  if (theCornerVertices[j] == theBoundaryVertices[i])
940  {
941  fixed_flags_boundary[i] = true;
942  num_fixed_boundary_flags++;
943  break;
944  }
945  }
946  }
947  printf("fixed %d of %d boundary vertices (those classified corner)\n",
948  num_fixed_boundary_flags, num_boundary_vertices);
949 
950 
951  // creates three intruction queues
952  Mesquite::InstructionQueue boundary_queue;
953  Mesquite::InstructionQueue interior_queue;
954 
955 
956  boundary_queue.set_slaved_ho_node_mode(Settings::SLAVE_ALL);
957  interior_queue.set_slaved_ho_node_mode(Settings::SLAVE_ALL);
958 
959  TShapeB1 targetMetric;
960 
961  IdealShapeTarget tc;
962 
963  TQualityMetric metric( &tc, &targetMetric );
964 
965  Mesquite::LPtoPTemplate* obj_func = new Mesquite::LPtoPTemplate(&metric, pOrder,
966  err);
967  if (MSQ_CHKERR(err)) {mfem::out << err << std::endl; exit(EXIT_FAILURE);}
968 
969  Mesquite::QuasiNewton* boundary_alg = new Mesquite::QuasiNewton( obj_func);
970  boundary_alg->use_element_on_vertex_patch();
971 
972  Mesquite::QuasiNewton* interior_alg = new Mesquite::QuasiNewton( obj_func );
973  interior_alg->use_global_patch();
974 
975  // **************Set stopping criterion**************
976  double grad_norm = 1e-5;
977  double successiveEps = 1e-5;
978  int boundary_outer = 5;
979  int boundary_inner = 3;
980 
981  // for boundary
982  Mesquite::TerminationCriterion* boundaryTermInner = new
983  Mesquite::TerminationCriterion();
984  Mesquite::TerminationCriterion* boundaryTermOuter = new
985  Mesquite::TerminationCriterion();
986  // boundaryTermInner->add_absolute_gradient_L2_norm(grad_norm);
987  // boundaryTermInner->add_relative_successive_improvement(successiveEps);
988  boundaryTermOuter->add_relative_successive_improvement(successiveEps);
989  boundaryTermOuter->add_iteration_limit(boundary_outer);
990  boundaryTermInner->add_iteration_limit(boundary_inner);
991 
992  ostringstream bndryStream;
993  bndryStream<<"boundary_"<<targetMetric.get_name()<<"_p"<<pOrder;
994  boundaryTermOuter->write_mesh_steps(bndryStream.str().c_str());
995  boundary_alg->set_outer_termination_criterion(boundaryTermOuter);
996  boundary_alg->set_inner_termination_criterion(boundaryTermInner);
997 
998  // for interior
999  Mesquite::TerminationCriterion* interiorTermInner = new
1000  Mesquite::TerminationCriterion();
1001  Mesquite::TerminationCriterion* interiorTermOuter = new
1002  Mesquite::TerminationCriterion();
1003  interiorTermInner->add_absolute_gradient_L2_norm(grad_norm);
1004  interiorTermInner->add_relative_successive_improvement(successiveEps);
1005  // interiorTermInner->add_iteration_limit(3); // for element_on_vertex_patch mode
1006  interiorTermInner->add_iteration_limit(100); // for global_patch mode
1007 
1008  ostringstream interiorStream;
1009  interiorStream<<"interior_"<<targetMetric.get_name()<<"_p"<<pOrder;
1010  interiorTermOuter->write_mesh_steps(interiorStream.str().c_str());
1011  interiorTermOuter->add_iteration_limit(1);
1012  interior_alg->set_outer_termination_criterion(interiorTermOuter);
1013  interior_alg->set_inner_termination_criterion(interiorTermInner);
1014 
1015  // ConditionNumberQualityMetric qm_metric;
1016  // QualityAssessor boundary_assessor,interior_assessor;
1017 
1018  // boundary_assessor.add_quality_assessment( &metric, 10 );
1019  // boundary_assessor.add_quality_assessment( &qm_metric );
1020 
1021  // interior_assessor.add_quality_assessment( &metric, 10 );
1022  // interior_assessor.add_quality_assessment( &qm_metric );
1023 
1024 
1025  // set the boundary instruction queue
1026  // boundary_queue.add_quality_assessor( &boundary_assessor, err );
1027  boundary_queue.set_master_quality_improver(boundary_alg, err);
1028  // boundary_queue.add_quality_assessor( &boundary_assessor, err );
1029  if (MSQ_CHKERR(err)) {mfem::out << err << std::endl; exit(EXIT_FAILURE);}
1030 
1031  // set the interior instruction queue
1032  // interior_queue.add_quality_assessor( &interior_assessor, err );
1033  interior_queue.set_master_quality_improver(interior_alg, err);
1034  // interior_queue.add_quality_assessor( &interior_assessor, err );
1035  if (MSQ_CHKERR(err)) {mfem::out << err << std::endl; exit(EXIT_FAILURE);}
1036 
1037  err.clear();
1038 
1039  std::vector<bool> fixed_flags(num_vertices);
1040 
1041  for (int j=0; j<mNumInterfaceSmoothIters; j++)
1042  {
1043 
1044  mfem::out<<" Boundary + Interior smoothing pass "<< j<<"....."<<endl;
1045 
1046  // smooth boundary only
1047  for (int i = 0; i < num_vertices; i++) { fixed_flags[i] = true; }
1048  mesh.vertices_set_fixed_flag(&(vertices[0]),fixed_flags,num_vertices, err);
1049  if (MSQ_CHKERR(err)) {mfem::out << err << std::endl; exit(EXIT_FAILURE);}
1050 
1051  mesh.vertices_set_fixed_flag(&(theBoundaryVertices[0]),fixed_flags_boundary,
1052  num_boundary_vertices, err);
1053  if (MSQ_CHKERR(err)) {mfem::out << err << std::endl; exit(EXIT_FAILURE);}
1054 
1055  // -----------------------------------------------------
1056  // debug
1057  //
1058  mesh.vertices_get_fixed_flag(&(vertices[0]),fixed_flags,num_vertices, err);
1059  if (MSQ_CHKERR(err)) {mfem::out << err << std::endl; exit(EXIT_FAILURE);}
1060  int num_fixed = 0;
1061  for (int i = 0; i < num_vertices; i++)
1062  {
1063  if (fixed_flags[i]) { num_fixed++; }
1064  }
1065 
1066  mfem::out << " For Boundary smooth, mesh has " << num_vertices <<
1067  " vertices and " << num_fixed << " are fixed. "<<endl;
1068  //
1069  // debug
1070  // -----------------------------------------------------
1071 
1072 
1073  boundary_queue.run_instructions(&mesh, mesh_domain, err);
1074  if (MSQ_CHKERR(err)) {mfem::out << err << std::endl; exit(EXIT_FAILURE);}
1075  mfem::out<<" boundary smooth completed in "<<boundaryTermOuter->get_iteration_count()
1076  <<" outer and "<<boundaryTermInner->get_iteration_count()
1077  <<" inner iterations."<<endl;
1078 
1079 
1080  // smooth interior only
1081 
1082  if (!mFixBndryInInterfaceSmooth)
1083  {
1084  //
1085  // let all interior vertices float during the interior smooth.
1086  //
1087  for (int i = 0; i < num_vertices; i++) { fixed_flags[i] = false; }
1088  mesh.vertices_set_fixed_flag(&(vertices[0]),fixed_flags,num_vertices, err);
1089  }
1090  else
1091  {
1092  //
1093  // use the app_fixed settings for the fixed state of boundary vertices.
1094  //
1095  mesh.vertices_set_fixed_flag(&(vertices[0]),app_fixed,num_vertices, err);
1096  }
1097  if (MSQ_CHKERR(err)) {mfem::out << err << std::endl; exit(EXIT_FAILURE);}
1098  for (int i = 0; i < num_boundary_vertices; i++) { fixed_flags[i] = true; }
1099 
1100  mesh.vertices_set_fixed_flag(&(theBoundaryVertices[0]),fixed_flags,
1101  num_boundary_vertices, err);
1102  if (MSQ_CHKERR(err)) {mfem::out << err << std::endl; exit(EXIT_FAILURE);}
1103 
1104  // -----------------------------------------------------
1105  // debug
1106  //
1107  mesh.vertices_get_fixed_flag(&(vertices[0]),fixed_flags,num_vertices, err);
1108  if (MSQ_CHKERR(err)) {mfem::out << err << std::endl; exit(EXIT_FAILURE);}
1109  num_fixed = 0;
1110  for (int i = 0; i < num_vertices; i++)
1111  {
1112  if (fixed_flags[i]) { num_fixed++; }
1113  }
1114 
1115  mfem::out << " For Interior smooth, mesh has " << num_vertices <<
1116  " vertices and " << num_fixed << " are fixed. "<<endl;
1117  //
1118  // debug
1119  // -----------------------------------------------------
1120 
1121  interior_queue.run_instructions(&mesh, &geom, err);
1122  if (MSQ_CHKERR(err)) {mfem::out << err << std::endl; exit(EXIT_FAILURE);}
1123  mfem::out<<" interior smooth completed in "<<interiorTermOuter->get_iteration_count()
1124  <<" outer and "<<interiorTermInner->get_iteration_count()
1125  <<" inner iterations."<<endl;
1126 
1127  }
1128 
1129 
1130  delete mesh_domain;
1131  delete interiorTermOuter;
1132  delete interiorTermInner;
1133  delete boundaryTermOuter;
1134  delete boundaryTermInner;
1135  delete interior_alg;
1136  delete boundary_alg;
1137  delete obj_func;
1138 
1139 }
1140 
1141 
1142 // Implementation of Mesh::MesquiteSmooth method
1143 
1144 void mfem::Mesh::MesquiteSmooth(const int mesquite_option)
1145 {
1146  mfem::MesquiteMesh msq_mesh(this);
1147  MsqDebug::enable(1);
1148  MsqPrintError err(mfem::err);
1149 
1150  Wrapper *method;
1151 
1152  const double vert_move_tol = 1e-3;
1153 
1154  switch (mesquite_option)
1155  {
1156  case 0: method = new LaplaceWrapper(); break;
1157  case 1: method = new UntangleWrapper(); break;
1158  case 2: method = new ShapeImprover(); break;
1159  case 3: method = new PaverMinEdgeLengthWrapper(vert_move_tol); break;
1160  }
1161 
1162  if ( mesquite_option < 4 )
1163  {
1164  // Specify SLAVE_NONE for high order node positions
1165  method->set_slaved_ho_node_mode(Settings::SLAVE_NONE);
1166 
1167  if ( this->Dimension() == 3 )
1168  {
1169  method->run_instructions(&msq_mesh, err);
1170  }
1171  else
1172  {
1173  Vector3D normal(0,0,1);
1174  Vector3D point(0,0,0);
1175  PlanarDomain mesh_plane(normal, point);
1176 
1177  method->run_instructions(&msq_mesh, &mesh_plane, err);
1178  }
1179  }
1180  else
1181  {
1182  // boundary perserving smoothing doesn't have a wrapper yet.
1183  BoundaryPreservingOptimization( msq_mesh );
1184  }
1185 }
1186 
1187 }
1188 
1189 #endif
void vertices_set_byte(const VertexHandle *vert_array, const unsigned char *byte_array, size_t array_size, MsqError &err)
Definition: mesquite.cpp:569
void get_all_vertices(std::vector< VertexHandle > &vertices, MsqError &err)
Definition: mesquite.cpp:516
void vertices_set_fixed_flag(const VertexHandle vert_array[], const std::vector< bool > &fixed_flag_array, size_t num_vtx, MsqError &err)
Definition: mesquite.cpp:613
TagHandle tag_create(const std::string &tag_name, TagType type, unsigned length, const void *default_value, MsqError &err)
Definition: mesquite.cpp:750
void vertex_set_coordinates(VertexHandle vertex, const Vector3D &coordinates, MsqError &err)
Definition: mesquite.cpp:550
int Size() const
Return the logical size of the array.
Definition: array.hpp:134
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:537
void BuildElementToDofTable() const
Definition: fespace.cpp:304
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:849
const Geometry::Type geom
Definition: ex1.cpp:40
Element::Type GetElementType(int i) const
Returns the type of element i.
Definition: mesh.cpp:5748
MesquiteMesh(mfem::Mesh *mfem_mesh)
Definition: mesquite.cpp:456
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:191
void tag_get_element_data(TagHandle handle, size_t num_elems, const ElementHandle *elem_array, void *tag_data, MsqError &err)
Definition: mesquite.cpp:807
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:846
void get_all_elements(std::vector< ElementHandle > &elements, MsqError &err)
Definition: mesquite.cpp:506
virtual void GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element &#39;elem&#39;.
Definition: fespace.cpp:2298
void vertices_get_coordinates(const VertexHandle vert_array[], MsqVertex *coordinates, size_t num_vtx, MsqError &err)
Definition: mesquite.cpp:526
void tag_properties(TagHandle handle, std::string &name_out, TagType &type_out, unsigned &length_out, MsqError &err)
Definition: mesquite.cpp:779
void vertex_set_byte(VertexHandle vertex, unsigned char byte, MsqError &err)
Definition: mesquite.cpp:561
void MesquiteSmooth(const int mesquite_option=0)
Definition: mesquite.cpp:1144
void vertices_get_attached_elements(const VertexHandle *vertex_array, size_t num_vertex, std::vector< ElementHandle > &elements, std::vector< size_t > &offsets, MsqError &err)
Definition: mesquite.cpp:661
void tag_get_vertex_data(TagHandle handle, size_t num_elems, const VertexHandle *node_array, void *tag_data, MsqError &err)
Definition: mesquite.cpp:833
void vertices_get_fixed_flag(const VertexHandle vert_array[], std::vector< bool > &fixed_flag_array, size_t num_vtx, MsqError &err)
Definition: mesquite.cpp:600
void tag_set_vertex_data(TagHandle handle, size_t num_elems, const VertexHandle *node_array, const void *tag_data, MsqError &err)
Definition: mesquite.cpp:820
void tag_delete(TagHandle handle, MsqError &err)
Definition: mesquite.cpp:763
const Element * GetElement(int i) const
Definition: mesh.hpp:942
virtual void GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element &#39;bel&#39;.
Definition: fespace.cpp:2418
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:629
int Dimension() const
Definition: mesh.hpp:911
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:71
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:414
void vertices_get_byte(const VertexHandle *vertex, unsigned char *byte_array, size_t array_size, MsqError &err)
Definition: mesquite.cpp:588
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:843
void GetNode(int i, double *coord) const
Definition: mesh.cpp:7292
int index(int i, int j, int nx, int ny)
Definition: life.cpp:237
void vertex_get_byte(const VertexHandle vertex, unsigned char *byte, MsqError &err)
Definition: mesquite.cpp:581
void elements_get_topologies(const ElementHandle *element_handle_array, EntityTopology *element_topologies, size_t num_elements, MsqError &err)
Definition: mesquite.cpp:684
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7343
void SetNode(int i, const double *coord)
Definition: mesh.cpp:7311
const Table & GetElementToDofTable() const
Return a reference to the internal Table that stores the lists of scalar dofs, for each mesh element...
Definition: fespace.hpp:700
int RowSize(int i) const
Definition: table.hpp:108
TagHandle tag_get(const std::string &name, MsqError &err)
Definition: mesquite.cpp:768
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
void tag_set_element_data(TagHandle handle, size_t num_elems, const ElementHandle *elem_array, const void *tag_data, MsqError &err)
Definition: mesquite.cpp:794
int GetAttribute(int i) const
Return the attribute of element i.
Definition: mesh.hpp:1197
void elements_get_attached_vertices(const ElementHandle *elem_handles, size_t num_elems, std::vector< VertexHandle > &vert_handles, std::vector< size_t > &offsets, MsqError &err)
Definition: mesquite.cpp:625
void GetBdrElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of boundary element i.
Definition: mesh.hpp:1015
int get_geometric_dimension(MsqError &err)
Definition: mesquite.cpp:501
Table * GetVertexToElementTable()
The returned Table must be destroyed by the caller.
Definition: mesh.cpp:5599