50 v(0) = p(0) + s*p(1) + s*p(2);
51 v(1) = p(1) + s*p(2) + s*p(0);
54 else if (p.
Size() == 2)
70 const double x = p(0), y = p(1);
72 return std::max(std::max(x - 0.25, -y), y - 1.0);
81 for (
int p = 0; p < np; p++)
84 fname << mesh_prefix <<
'.' << setfill(
'0') << setw(6) << p;
88 cerr <<
"Can not open mesh file: " << fname.str().c_str()
90 for (p--; p >= 0; p--)
96 mesh_array[p] =
new Mesh(meshin, 1, 0);
100 for (
int i = 0; i < mesh_array[p]->GetNE(); i++)
102 mesh_array[p]->GetElement(i)->SetAttribute(p+1);
104 for (
int i = 0; i < mesh_array[p]->GetNBE(); i++)
106 mesh_array[p]->GetBdrElement(i)->SetAttribute(p+1);
110 mesh =
new Mesh(mesh_array, np);
112 for (
int p = 0; p < np; p++)
114 delete mesh_array[np-1-p];
127 for (
int i = 0; i < mesh->
GetNBE(); i++)
132 for (
int j = 0; j < nv; j++)
138 for (
int i = 0; i < v2v.
Size(); i++)
152 for (
int i = 0; i < v2v.
Size(); i++)
164 for (
int i = 0; i < mesh->
GetNBE(); i++)
170 for (
int j = 0; j < nv; j++)
198 if (dynamic_cast<const H1_FECollection*>(fec))
212 for (
int i=0; i<mesh->
GetNBE(); i++)
215 mesh->
GetNodes()->GetSubVector(vdofs, v);
217 fes_copy->GetElementVDofs(i, bvdofs);
223 cout <<
"\nDiscontinuous nodes not yet supported" << endl;
230 int main (
int argc,
char *argv[])
233 const char *mesh_file =
"../../data/beam-hex.mesh";
237 args.
AddOption(&mesh_file,
"-m",
"--mesh",
238 "Mesh file to visualize.");
240 "Load mesh from multiple processors.");
241 args.
AddOption(&refine,
"-ref",
"--refinement",
"-no-ref",
"--no-refinement",
242 "Prepare the mesh for refinement or not.");
251 cout <<
"Visualize and manipulate a serial mesh:\n"
252 <<
" mesh-explorer -m <mesh_file>\n"
253 <<
"Visualize and manipulate a parallel mesh:\n"
254 <<
" mesh-explorer -np <#proc> -m <mesh_prefix>\n" << endl
262 Mesh *bdr_mesh = NULL;
265 mesh =
new Mesh(mesh_file, 1, refine);
297 cout <<
"boundary attribs :";
302 cout <<
'\n' <<
"material attribs :";
308 cout <<
"mesh curvature : ";
315 cout <<
"NONE" << endl;
320 cout <<
"What would you like to do?\n"
322 "c) Change curvature\n"
327 "m) View materials\n"
330 "h) View element sizes, h\n"
331 "k) View element ratios, kappa\n"
332 "x) Print sub-element stats\n"
333 "f) Find physical point in reference space\n"
334 "p) Generate a partitioning\n"
335 "o) Reorder elements\n"
336 "S) Save in MFEM format\n"
337 "V) Save in VTK format (only linear and quadratic meshes)\n"
340 "Z) Save in MFEM format with compression\n"
355 "Choose type of refinement:\n"
356 "s) standard refinement with Mesh::UniformRefinement()\n"
357 "b) Mesh::UniformRefinement() (bisection for tet meshes)\n"
358 "u) uniform refinement with a factor\n"
359 "g) non-uniform refinement (Gauss-Lobatto) with a factor\n"
360 "l) refine locally using the region() function\n"
377 cout <<
"enter refinement factor --> " << flush;
380 if (ref_factor <= 1 || ref_factor > 32) {
break; }
383 Mesh *rmesh =
new Mesh(mesh, ref_factor, ref_type);
392 for (
int i = 0; i < mesh->
GetNE(); i++)
402 marked_elements.
Append(i);
417 cout <<
"enter new order for mesh curvature --> " << flush;
426 cout <<
"scaling factor ---> " << flush;
432 for (
int i = 0; i < mesh->
GetNV(); i++)
460 cout <<
"jitter factor ---> " << flush;
467 cerr <<
"The mesh should have nodes, introduce curvature first!\n";
483 for (
int i = 0; i < fespace->
GetNE(); i++)
486 for (
int j = 0; j < dofs.
Size(); j++)
494 for (
int i = 0; i < fespace->
GetNDofs(); i++)
496 for (
int d = 0; d <
dim; d++)
503 cout <<
"move boundary nodes? [y/n] ---> " << flush;
510 for (
int i = 0; i < fespace->
GetNBE(); i++)
513 for (
int j = 0; j < vdofs.
Size(); j++)
530 double min_det_J, max_det_J, min_det_J_z, max_det_J_z;
531 double min_kappa, max_kappa, max_ratio_det_J_z;
533 max_det_J = max_kappa = max_ratio_det_J_z = -
infinity();
534 cout <<
"subdivision factor ---> " << flush;
537 bad_elems_by_geom = 0;
538 for (
int i = 0; i < mesh->
GetNE(); i++)
553 double det_J = J.
Det();
557 min_det_J_z = fmin(min_det_J_z, det_J);
558 max_det_J_z = fmax(max_det_J_z, det_J);
560 min_kappa = fmin(min_kappa, kappa);
561 max_kappa = fmax(max_kappa, kappa);
564 fmax(max_ratio_det_J_z, max_det_J_z/min_det_J_z);
565 min_det_J = fmin(min_det_J, min_det_J_z);
566 max_det_J = fmax(max_det_J, max_det_J_z);
567 if (min_det_J_z <= 0.0)
570 bad_elems_by_geom[
geom]++;
573 cout <<
"\nbad elements = " << nz;
579 cout <<
"\nmin det(J) = " << min_det_J
580 <<
"\nmax det(J) = " << max_det_J
581 <<
"\nglobal ratio = " << max_det_J/min_det_J
582 <<
"\nmax el ratio = " << max_ratio_det_J_z
583 <<
"\nmin kappa = " << min_kappa
584 <<
"\nmax kappa = " << max_kappa << endl;
590 cout <<
"\npoint in physical space ---> " << flush;
591 for (
int i = 0; i < sdim; i++)
593 cin >> point_mat(i,0);
601 cout <<
"point in reference space:";
602 if (elem_ids[0] == -1)
604 cout <<
" NOT FOUND!\n";
608 cout <<
" element " << elem_ids[0] <<
", ip =";
609 cout <<
" " << ips[0].x;
612 cout <<
" " << ips[0].y;
615 cout <<
" " << ips[0].z;
624 cout <<
"What type of reordering?\n"
625 "h) Hilbert spatial sort\n"
640 if (mk ==
'm' || mk ==
'b' || mk ==
'e' || mk ==
'v' || mk ==
'h' ||
641 mk ==
'k' || mk ==
'p')
653 for (
int i = 0; i < mesh->
GetNE(); i++)
668 bdr_attr.
SetSpace(bdr_attr_fespace);
669 for (
int i = 0; i < bdr_mesh->
GetNE(); i++)
671 bdr_part[i] = (bdr_attr(i) = bdr_mesh->
GetAttribute(i)) - 1;
689 double a = double(rand()) / (double(RAND_MAX) + 1.);
690 int el0 = (int)floor(a * mesh->
GetNE());
691 cout <<
"Generating coloring starting with element " << el0+1
692 <<
" / " << mesh->
GetNE() << endl;
694 for (
int i = 0; i < coloring.
Size(); i++)
696 attr(i) = coloring[i];
698 cout <<
"Number of colors: " << attr.
Max() + 1 << endl;
699 for (
int i = 0; i < mesh->
GetNE(); i++)
702 attr(i) = part[i] = i;
712 for (
int i = 0; i < mesh->
GetNE(); i++)
722 attr(i) = -pow(-attr(i), 1.0/
double(dim));
726 attr(i) = pow(attr(i), 1.0/
double(dim));
728 h_min = min(h_min, attr(i));
729 h_max = max(h_max, attr(i));
731 cout <<
"h_min = " << h_min <<
", h_max = " << h_max << endl;
737 for (
int i = 0; i < mesh->
GetNE(); i++)
749 int *partitioning = NULL, np;
750 cout <<
"What type of partitioning?\n"
752 "s) Simple 1D split of the element sequence\n"
753 "0) METIS_PartGraphRecursive (sorted neighbor lists)\n"
754 "1) METIS_PartGraphKway (sorted neighbor lists)"
756 "2) METIS_PartGraphVKway (sorted neighbor lists)\n"
757 "3) METIS_PartGraphRecursive\n"
758 "4) METIS_PartGraphKway\n"
759 "5) METIS_PartGraphVKway\n"
766 cout <<
"Enter nx: " << flush;
767 cin >> nxyz[0]; np = nxyz[0];
770 cout <<
"Enter ny: " << flush;
771 cin >> nxyz[1]; np *= nxyz[1];
774 cout <<
"Enter nz: " << flush;
775 cin >> nxyz[2]; np *= nxyz[2];
782 cout <<
"Enter number of processors: " << flush;
785 partitioning =
new int[mesh->
GetNE()];
786 for (
int i = 0; i < mesh->
GetNE(); i++)
788 partitioning[i] = i * np / mesh->
GetNE();
793 int part_method = pk -
'0';
794 if (part_method < 0 || part_method > 5)
798 cout <<
"Enter number of processors: " << flush;
804 const char part_file[] =
"partitioning.txt";
805 ofstream opart(part_file);
806 opart <<
"number_of_elements " << mesh->
GetNE() <<
'\n'
807 <<
"number_of_processors " << np <<
'\n';
808 for (
int i = 0; i < mesh->
GetNE(); i++)
810 opart << partitioning[i] <<
'\n';
812 cout <<
"Partitioning file: " << part_file << endl;
816 for (
int i = 0; i < mesh->
GetNE(); i++)
818 proc_el[partitioning[i]]++;
820 int min_el = proc_el[0], max_el = proc_el[0];
821 for (
int i = 1; i < np; i++)
823 if (min_el > proc_el[i])
827 if (max_el < proc_el[i])
832 cout <<
"Partitioning stats:\n"
834 << setw(12) <<
"minimum"
835 << setw(12) <<
"average"
836 << setw(12) <<
"maximum"
837 << setw(12) <<
"total" <<
'\n';
839 << setw(12) << min_el
840 << setw(12) << double(mesh->
GetNE())/np
841 << setw(12) << max_el
842 << setw(12) << mesh->
GetNE() << endl;
849 for (
int i = 0; i < mesh->
GetNE(); i++)
851 attr(i) = part[i] = partitioning[i];
853 delete [] partitioning;
856 char vishost[] =
"localhost";
861 sol_sock.precision(14);
864 sol_sock <<
"fem2d_gf_data_keys\n";
867 mesh->
Print(sol_sock);
874 mesh->
Print(sol_sock);
875 for (
int i = 0; i < mesh->
GetNE(); i++)
886 sol_sock <<
"RjlmAb***********";
898 sol_sock <<
"fem3d_gf_data_keys\n";
899 if (mk ==
'v' || mk ==
'h' || mk ==
'k')
901 mesh->
Print(sol_sock);
905 bdr_mesh->
Print(sol_sock);
906 bdr_attr.
Save(sol_sock);
909 sol_sock <<
"pppppp" <<
"pppppp" <<
"pppppp";
916 mesh->
Print(sol_sock);
917 for (
int i = 0; i < mesh->
GetNE(); i++)
945 cout <<
"Unable to connect to "
946 << vishost <<
':' << visport << endl;
949 delete bdr_attr_fespace;
954 const char mesh_file[] =
"mesh-explorer.mesh";
955 ofstream omesh(mesh_file);
958 cout <<
"New mesh file: " << mesh_file << endl;
963 const char mesh_file[] =
"mesh-explorer.vtk";
964 ofstream omesh(mesh_file);
967 cout <<
"New VTK mesh file: " << mesh_file << endl;
973 const char mesh_file[] =
"mesh-explorer.mesh.gz";
977 cout <<
"New mesh file: " << mesh_file << endl;
int GetNPoints() const
Returns the number of the points in the integration rule.
Geometry::Type GetGeometryType() const
int Size() const
Logical size of the array.
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
virtual void Print(std::ostream &out=mfem::out) const
int * CartesianPartitioning(int nxyz[])
int GetNDofs() const
Returns number of degrees of freedom.
Class for an integration rule - an Array of IntegrationPoint.
const double * GetVertex(int i) const
Return pointer to vertex i's coordinates.
Class for grid function - Vector with associated FE space.
int DofToVDof(int dof, int vd, int ndofs=-1) const
Mesh * read_par_mesh(int np, const char *mesh_prefix)
virtual void GetVertices(Array< int > &v) const =0
Returns element's vertices.
static void PrintElementsByGeometry(int dim, const Array< int > &num_elems_by_geom, std::ostream &out)
Auxiliary method used by PrintCharacteristics().
void JacToPerfJac(int GeomType, const DenseMatrix &J, DenseMatrix &PJ) const
int * GeneratePartitioning(int nparts, int part_method=1)
int GetNBE() const
Returns number of boundary elements.
void PrintWithPartitioning(int *partitioning, std::ostream &out, int elem_attr=0) const
Prints the mesh with boundary elements given by the boundary of the subdomains, so that the boundary ...
void MakeOwner(FiniteElementCollection *_fec)
Make the GridFunction the owner of fec and fes.
void NewNodes(GridFunction &nodes, bool make_owner=false)
Replace the internal node GridFunction with the given GridFunction.
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Data type dense matrix using column-major storage.
int Size() const
Returns the size of the vector.
void Transform(void(*f)(const Vector &, Vector &))
int GetNE() const
Returns number of elements.
void PrintHelp(std::ostream &out) const
void Randomize(int seed=0)
Set random values in the vector.
const IntegrationPoint & GetCenter(int GeomType)
Return the center of the given Geometry::Type, GeomType.
int main(int argc, char *argv[])
Geometry::Type GetElementBaseGeometry(int i) const
void DeleteAll()
Delete whole array.
Mesh * skin_mesh(Mesh *mesh)
void AddVertex(const double *)
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
int GetNE() const
Returns number of elements in the mesh.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
int GetNBE() const
Returns number of boundary elements in the mesh.
void AddSegment(const int *vi, int attr=1)
int Append(const T &el)
Append element to array, resize if necessary.
void transformation(const Vector &p, Vector &v)
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
virtual void GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
GeometryRefiner GlobGeometryRefiner
int GetAttribute() const
Return element's attribute.
void GetHilbertElementOrdering(Array< int > &ordering)
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Nodes: x_i = i/(n-1), i=0,...,n-1.
FiniteElementSpace * FESpace()
double GetElementSize(int i, int type=0)
Get the size of the i-th element relative to the perfect reference element.
int SpaceDimension() const
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
void PrintVTK(std::ostream &out)
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
void AddQuad(const int *vi, int attr=1)
void GetColumnReference(int c, Vector &col)
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name...
NURBSExtension * NURBSext
Optional NURBS mesh extension.
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
void PrintError(std::ostream &out) const
virtual const char * Name() const
void AddTriangle(const int *vi, int attr=1)
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
double Max() const
Returns the maximal element of the vector.
const FiniteElementSpace * GetNodalFESpace() const
virtual void Finalize(bool refine=false, bool fix_orientation=false)
Finalize the construction of a general Mesh.
void ReorderElements(const Array< int > &ordering, bool reorder_vertices=true)
void PrintOptions(std::ostream &out) const
double CalcSingularvalue(const int i) const
Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
virtual int GetNVertices() const =0
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
virtual int FindPoints(DenseMatrix &point_mat, Array< int > &elem_ids, Array< IntegrationPoint > &ips, bool warn=true, InverseElementTransformation *inv_trans=NULL)
Find the ids of the elements that contain the given points, and their corresponding reference coordin...
const FiniteElementCollection * FEColl() const
void GetNodes(Vector &node_coord) const
void GetElementColoring(Array< int > &colors, int el0=0)
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i'th boundary element.
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
double region(const Vector &p)
Abstract data type element.
int GetAttribute(int i) const
Return the attribute of element i.
Array< int > attributes
A list of all unique element attributes used by the Mesh.
const Element * GetBdrElement(int i) const
void PrintCharacteristics(Vector *Vh=NULL, Vector *Vk=NULL, std::ostream &out=mfem::out)
Compute and print mesh characteristics such as number of vertices, number of elements, number of boundary elements, minimal and maximal element sizes, minimal and maximal element aspect ratios, etc.