static char help[] = "Tests the DMShellCreate() and PetscSectionCreate() for fieldsplit initialization.";

#include "petscsys.h"
#include "petscvec.h"
#include "petscmat.h"
#include "petscdm.h"
#include "petscdmshell.h"
#include "petscksp.h"

int main(int argc, char **argv)
{
  int               i, indsfield[8];
  PetscInt          numeq = 8, start = 1, NF, NE, numFields = 2, inds[8];
  PetscErrorCode    ierr;
  DM                FSDM;
  PetscSection      FSSection;
  Mat               Kmat;
  Vec               rhs, sol;
  KSP               kspsol;
  const PetscScalar K[] = {4,2,0,0,2,1,0,0,2,8,2,0,1,4,0,0,0,2,8,2,0,0,3,2,0,0,2,4,0,0,2,1,0,0,0,0,1,0,0,0,0,0,0,0,0,3,0,0,0,0,0,0,0,0,2,0,0,0,0,0,0,0,0,1};
  const PetscScalar r[] = {3,3,3,3,3,3,3,3};
  
  indsfield[0]  = 0;
  indsfield[1]  = 0;
  indsfield[2]  = 0;
  indsfield[3]  = 0;
  indsfield[4]  = 1;
  indsfield[5]  = 1;
  indsfield[6]  = 1;
  indsfield[7]  = 1;
  
  for(i=0; i<8; i++)
    {
      inds[i] = i;
    }
  
  // We would like to use a DMShell to represent this data structure.
  // To do so we start by creating a section with 10points having one dof each.
  ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);
  
  ierr = PetscSectionCreate(PETSC_COMM_WORLD, &FSSection);CHKERRQ(ierr);
  ierr = PetscSectionSetNumFields(FSSection, numFields);CHKERRQ(ierr);  // <---- We create the two fields
  ierr = PetscSectionSetChart(FSSection, start, numeq+1);CHKERRQ(ierr); // <---- We set the chart for the section /!\ This needs to be done after the creation of the fields
  ierr = PetscSectionGetChart(FSSection, &NF, &NE);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD, "start=%i, End=%i.\n",NF,NE);CHKERRQ(ierr);
  
  for(i=0; i<numeq; i++)
    {
      ierr = PetscSectionSetDof(FSSection, i+1, 1);CHKERRQ(ierr);
      ierr = PetscSectionSetFieldDof(FSSection, i+1, indsfield[i], 1);CHKERRQ(ierr);  // <---- We assign the dofs to their field
    }
  ierr = PetscSectionSetUp(FSSection);CHKERRQ(ierr);
  ierr = PetscSectionView(FSSection, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  
  
  ierr = DMShellCreate(PETSC_COMM_WORLD, &FSDM);CHKERRQ(ierr);
  ierr = DMSetDefaultSection(FSDM, FSSection);CHKERRQ(ierr);
  ierr = DMSetUp(FSDM);CHKERRQ(ierr);
  //  ierr = DMView(FSDM, PETSC_VIEWER_STDOUT_WORLD);
  
  /* Create and Assemble the matrix K */
  ierr = MatCreate(PETSC_COMM_WORLD, &Kmat);CHKERRQ(ierr);
  ierr = MatSetSizes(Kmat, numeq, numeq, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr);
  ierr = MatSetType(Kmat, MATAIJ);CHKERRQ(ierr);
  ierr = MatSetFromOptions(Kmat);CHKERRQ(ierr);
  ierr = MatSetUp(Kmat);CHKERRQ(ierr);
  if(ierr != 0)
    {
      printf("There has been a problem creating Kmat.\n");
      return 1;
    }
  
  ierr = MatZeroEntries(Kmat);CHKERRQ(ierr);
  ierr = MatSetValues(Kmat, 8, inds, 8, inds, K, ADD_VALUES);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(Kmat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(Kmat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  if(ierr != 0)
    {
      printf("There has been a problem assembling Kmat.\n");
      return 1;
    }
  /* ierr = MatView(Kmat, PETSC_VIEWER_STDOUT_WORLD); */
  
  /* Create and Assemble the right hand side and the solution vector */
  ierr = VecCreate(PETSC_COMM_WORLD, &rhs);CHKERRQ(ierr);
  ierr = VecSetSizes(rhs, numeq, PETSC_DECIDE);CHKERRQ(ierr);
  ierr = VecSetFromOptions(rhs);CHKERRQ(ierr);
  ierr = VecSetUp(rhs);CHKERRQ(ierr);
  ierr = VecDuplicate(rhs, &sol);CHKERRQ(ierr);
  if(ierr != 0)
    {
      printf("There has been a problem creating rhs and sol.\n");
      return 1;
    }
  
  ierr = VecSetValues(rhs, 8, inds, r, ADD_VALUES);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(rhs);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(rhs);CHKERRQ(ierr);
  if(ierr != 0)
    {
      printf("There has been a problem assembling rhs.\n");
      return 1;
    }
  /* ierr = VecView(rhs, PETSC_VIEWER_STDOUT_WORLD); */
  
  /* Create the kspsol and pc */
  ierr = KSPCreate(PETSC_COMM_WORLD, &kspsol);CHKERRQ(ierr);
  ierr = KSPSetDM(kspsol, FSDM);CHKERRQ(ierr);
  ierr = KSPSetDMActive(kspsol, PETSC_FALSE);CHKERRQ(ierr);
  ierr = KSPSetOperators(kspsol, Kmat, Kmat, DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(kspsol);CHKERRQ(ierr);

  ierr = KSPSolve(kspsol, rhs, sol);CHKERRQ(ierr);
  printf("The solution to this problem is:\n");
  ierr = VecView(sol, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  
  ierr = KSPDestroy(&kspsol);CHKERRQ(ierr);
  ierr = VecDestroy(&sol);CHKERRQ(ierr);
  ierr = VecDestroy(&rhs);CHKERRQ(ierr);
  ierr = MatDestroy(&Kmat);CHKERRQ(ierr);
  ierr = DMDestroy(&FSDM);CHKERRQ(ierr);
  ierr = PetscSectionDestroy(&FSSection);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return 0;
}