1c4762a1bSJed Brown static char help[] = "Define a simple field over the mesh\n\n"; 2c4762a1bSJed Brown 3c4762a1bSJed Brown #include <petscdmplex.h> 4c4762a1bSJed Brown 5d71ae5a4SJacob Faibussowitsch int main(int argc, char **argv) 6d71ae5a4SJacob Faibussowitsch { 730602db0SMatthew G. Knepley DM dm; 8c4762a1bSJed Brown Vec u; 9c4762a1bSJed Brown PetscSection section; 10c4762a1bSJed Brown PetscViewer viewer; 1130602db0SMatthew G. Knepley PetscInt dim, numFields, numBC, i; 12c4762a1bSJed Brown PetscInt numComp[3]; 13c4762a1bSJed Brown PetscInt numDof[12]; 14c4762a1bSJed Brown PetscInt bcField[1]; 15c4762a1bSJed Brown IS bcPointIS[1]; 16c4762a1bSJed Brown 17327415f7SBarry Smith PetscFunctionBeginUser; 189566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 19c4762a1bSJed Brown /* Create a mesh */ 209566063dSJacob Faibussowitsch PetscCall(DMCreate(PETSC_COMM_WORLD, &dm)); 219566063dSJacob Faibussowitsch PetscCall(DMSetType(dm, DMPLEX)); 229566063dSJacob Faibussowitsch PetscCall(DMSetFromOptions(dm)); 239566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 249566063dSJacob Faibussowitsch PetscCall(DMGetDimension(dm, &dim)); 25c4762a1bSJed Brown /* Create a scalar field u, a vector field v, and a surface vector field w */ 26c4762a1bSJed Brown numFields = 3; 27c4762a1bSJed Brown numComp[0] = 1; 28c4762a1bSJed Brown numComp[1] = dim; 29c4762a1bSJed Brown numComp[2] = dim - 1; 30c4762a1bSJed Brown for (i = 0; i < numFields * (dim + 1); ++i) numDof[i] = 0; 31c4762a1bSJed Brown /* Let u be defined on vertices */ 32c4762a1bSJed Brown numDof[0 * (dim + 1) + 0] = 1; 33c4762a1bSJed Brown /* Let v be defined on cells */ 34c4762a1bSJed Brown numDof[1 * (dim + 1) + dim] = dim; 35c4762a1bSJed Brown /* Let w be defined on faces */ 36c4762a1bSJed Brown numDof[2 * (dim + 1) + dim - 1] = dim - 1; 37c4762a1bSJed Brown /* Setup boundary conditions */ 38c4762a1bSJed Brown numBC = 1; 39c4762a1bSJed Brown /* Prescribe a Dirichlet condition on u on the boundary 40c4762a1bSJed Brown Label "marker" is made by the mesh creation routine */ 41c4762a1bSJed Brown bcField[0] = 0; 429566063dSJacob Faibussowitsch PetscCall(DMGetStratumIS(dm, "marker", 1, &bcPointIS[0])); 43c4762a1bSJed Brown /* Create a PetscSection with this data layout */ 449566063dSJacob Faibussowitsch PetscCall(DMSetNumFields(dm, numFields)); 459566063dSJacob Faibussowitsch PetscCall(DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, §ion)); 469566063dSJacob Faibussowitsch PetscCall(ISDestroy(&bcPointIS[0])); 47c4762a1bSJed Brown /* Name the Field variables */ 489566063dSJacob Faibussowitsch PetscCall(PetscSectionSetFieldName(section, 0, "u")); 499566063dSJacob Faibussowitsch PetscCall(PetscSectionSetFieldName(section, 1, "v")); 509566063dSJacob Faibussowitsch PetscCall(PetscSectionSetFieldName(section, 2, "w")); 51*0cae8135SMatthew G. Knepley PetscCall(PetscSectionViewFromOptions(section, NULL, "-mysection_view")); 52c4762a1bSJed Brown /* Tell the DM to use this data layout */ 539566063dSJacob Faibussowitsch PetscCall(DMSetLocalSection(dm, section)); 54c4762a1bSJed Brown /* Create a Vec with this layout and view it */ 559566063dSJacob Faibussowitsch PetscCall(DMGetGlobalVector(dm, &u)); 569566063dSJacob Faibussowitsch PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer)); 579566063dSJacob Faibussowitsch PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK)); 589566063dSJacob Faibussowitsch PetscCall(PetscViewerFileSetName(viewer, "sol.vtu")); 599566063dSJacob Faibussowitsch PetscCall(VecView(u, viewer)); 609566063dSJacob Faibussowitsch PetscCall(PetscViewerDestroy(&viewer)); 619566063dSJacob Faibussowitsch PetscCall(DMRestoreGlobalVector(dm, &u)); 62c4762a1bSJed Brown /* Cleanup */ 639566063dSJacob Faibussowitsch PetscCall(PetscSectionDestroy(§ion)); 649566063dSJacob Faibussowitsch PetscCall(DMDestroy(&dm)); 659566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 66b122ec5aSJacob Faibussowitsch return 0; 67c4762a1bSJed Brown } 68c4762a1bSJed Brown 69c4762a1bSJed Brown /*TEST 70c4762a1bSJed Brown 71c4762a1bSJed Brown test: 72c4762a1bSJed Brown suffix: 0 73c4762a1bSJed Brown requires: triangle 74*0cae8135SMatthew G. Knepley args: -mysection_view -info :~sys,mat 75*0cae8135SMatthew G. Knepley 76c4762a1bSJed Brown test: 77c4762a1bSJed Brown suffix: 1 78c4762a1bSJed Brown requires: ctetgen 79*0cae8135SMatthew G. Knepley args: -dm_plex_dim 3 -mysection_view -info :~sys,mat 80*0cae8135SMatthew G. Knepley 81*0cae8135SMatthew G. Knepley test: 82*0cae8135SMatthew G. Knepley suffix: filter_0 83*0cae8135SMatthew G. Knepley args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3 -dm_plex_label_lower 0,1,2,3,4,6 \ 84*0cae8135SMatthew G. Knepley -dm_refine 1 -dm_plex_transform_type transform_filter -dm_plex_transform_active lower \ 85*0cae8135SMatthew G. Knepley -dm_view 86*0cae8135SMatthew G. Knepley 87*0cae8135SMatthew G. Knepley test: 88*0cae8135SMatthew G. Knepley suffix: submesh_0 89*0cae8135SMatthew G. Knepley requires: triangle 90*0cae8135SMatthew G. Knepley args: -dm_plex_label_middle 9,12,15,23,31 -dm_plex_submesh middle -dm_view 91c4762a1bSJed Brown 92c4762a1bSJed Brown TEST*/ 93