regrefine.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2010 Werner Scholz

    www:   http://www.magpar.net/
    email: magpar(at)magpar.net

    magpar is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    magpar is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with magpar; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

static char const Id[] = "$Id: regrefine.c 3037 2010-03-30 07:10:51Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "init/init.h"
#include "util/util.h"

int RegularRefinement(GridData *gdata)
{
  const int midnod[4][4]=
            {
              {C_UNK,4,6,7},
              {C_UNK,C_UNK,5,8},
              {C_UNK,C_UNK,C_UNK,9},
              {C_UNK,C_UNK,C_UNK,C_UNK}
            };
  const int child[8][4]=
            {
              {0,4,6,7},
              {4,1,5,8},
              {6,5,2,9},
              {7,8,9,3},
              {4,6,7,8},
              {4,5,6,8},
              {6,8,5,9},
              {6,7,8,9}
            };


  MagparFunctionLogBegin;

  int rank,size;
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

  int refine;
  PetscTruth flg;
  ierr = PetscOptionsGetInt(PETSC_NULL,"-refine",(PetscInt*)&refine,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    refine=0;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -refine not found, using default value: %i\n",
      refine
    );
  }

  PetscPrintf(PETSC_COMM_WORLD,
    "Requested number of refinement steps: %i\n",
    refine
  );

  if (rank) {
    for (int r=0;r<refine;r++) {
      Vec Mnew;
      ierr = VecCreate(PETSC_COMM_WORLD,&Mnew);CHKERRQ(ierr);
      ierr = VecSetSizes(Mnew,ND*gdata->ln_vert,PETSC_DECIDE);CHKERRQ(ierr);
      ierr = VecSetFromOptions(Mnew);CHKERRQ(ierr);
      ierr = VecSetBlockSize(Mnew,ND);CHKERRQ(ierr);

      ierr = VecDestroy(gdata->M);CHKERRQ(ierr);
      gdata->M=Mnew;

      PetscReal devNorm;
      devNorm=RenormVec(gdata->M,1.0,0.0,ND);
    }
    MagparFunctionLogReturn(0);
  }

  /* get connectivity matrix of nodes */
  /* recommended size (METIS manual):
     sizeof(ia): (n_vert+1)
     sizeof(ja): 15*n_vert
  */

  for (int r=0;r<refine;r++) {
    PetscPrintf(PETSC_COMM_WORLD,"Starting refinement step: %i\n",r+1);

    int *nxadj,*nadjncy;
    ierr = Mesh2Nodal(gdata->n_ele,gdata->n_vert,gdata->elevert,&nxadj,&nadjncy);CHKERRQ(ierr);

    /* the last entry in nxadj (nxadj[gdata->n_vert])
       gives us the total number of vertex neighbours,
       each pair of neighbours appears twice,
       thus we have nxadj[gdata->n_vert]/2 pairs of neighbours (=edges),
       thus the total number of edges is
       nxadj[gdata->n_vert]/2
    */

    int n_edge;
    n_edge=nxadj[gdata->n_vert]/2;
    PetscPrintf(PETSC_COMM_WORLD,
      "Number of edges: %i\n",
      n_edge
    );

    PetscPrintf(PETSC_COMM_WORLD,
      "Generating a FE mesh with %i nodes and %i elements.\n",
      gdata->n_vert+n_edge,
      8*gdata->n_ele
    );


    /* allocate memory for ids of new vertices */
    int *nvert;
    ierr = PetscMalloc(2*n_edge*sizeof(int),&nvert);CHKERRQ(ierr);

    for (int i=0;i<2*n_edge;i++) nvert[i]=C_UNK;

    /* reset counter for new vertices */
    int cnt;
    cnt=gdata->n_vert;

    /* allocate memory for cartesian coordinates of old and new vertices */
    PetscReal *nvertxyz;
    ierr = PetscMalloc(ND*(gdata->n_vert+n_edge)*sizeof(PetscReal),&nvertxyz);CHKERRQ(ierr);

    ierr = PetscMemcpy(nvertxyz,gdata->vertxyz,ND*gdata->n_vert*sizeof(PetscReal));CHKERRQ(ierr);

    /* allocate memory for new elements */
    int *nelevert;
    ierr = PetscMalloc(NV*8*gdata->n_ele*sizeof(int),&nelevert);CHKERRQ(ierr);

    /* regular refinement
       cf.
       J. Bey, Tetrahedral Grid Refinement. Computing 55 (1995) 355.
       http://www.igpm.rwth-aachen.de/bey/
       http://www.igpm.rwth-aachen.de/Download/reports/bey/tetra.ps.gz
       http://www.igpm.rwth-aachen.de/Download/reports/bey/agm.ps.gz
       p. 62
    */
    for (int i=0;i<gdata->n_ele;i++) {
      /* get vertex ids */
      int t_v[10];
      ierr = PetscMemcpy(t_v,gdata->elevert+i*NV,NV*sizeof(int));CHKERRQ(ierr);

      /* find new vertices */
      for (int m=0;m<NV-1;m++) {
        for (int n=m+1;n<NV;n++) {
          /* we take the minimum of the "parent node ids" since
             they are "equal partners" and have only one common child,
             and we want to have a unique address of the child's id in nvert
          */
          int lid,hid;
          lid=PetscMin(t_v[m],t_v[n]);
          hid=PetscMax(t_v[m],t_v[n]);

          int nnb;
          nnb=nxadj[lid+1]-nxadj[lid];

          /* search in list of newly created nodes first */

          /* go through the list of neigbours */
          int j;
          for (j=0;j<nnb;j++) {
            if (nadjncy[nxadj[lid]+j]==hid) {
              /* we found the neighbouring node */

              /* check if it the midnode has already been "created"
              */
              if (nvert[nxadj[lid]+j] == C_UNK) {
                nvert[nxadj[lid]+j]=cnt;
                /* set cartesian coordinates */
                for (int k=0;k<ND;k++) {
                  nvertxyz[ND*cnt+k]=
                    0.5*(gdata->vertxyz[ND*lid+k]+gdata->vertxyz[ND*hid+k]);
                }
                /* increase vertex counter */
                cnt++;
              }
              break;
            }
          }
          /* the neighbour must have been found
             (and the for loop terminated by the break command)
          */
          assert(j<nxadj[lid+1]-nxadj[lid]);

          t_v[midnod[m][n]]=nvert[nxadj[lid]+j];
        }
      }
      /* create new elements */
      for (int j=0;j<8;j++) {
        for (int k=0;k<NV;k++) {
          nelevert[NV*8*i+NV*j+k]=t_v[child[j][k]];
        }
      }
    }

    assert(cnt==gdata->n_vert+n_edge);

    ierr = PetscFree(gdata->vertxyz);CHKERRQ(ierr);
    gdata->vertxyz=nvertxyz;

    ierr = PetscFree(gdata->elevert);CHKERRQ(ierr);
    gdata->elevert=nelevert;

    /* interpolate magnetization on new vertices */

    Vec Mnew;
    ierr = VecCreate(PETSC_COMM_WORLD,&Mnew);CHKERRQ(ierr);
    ierr = VecSetSizes(Mnew,ND*(gdata->n_vert+n_edge),PETSC_DECIDE);CHKERRQ(ierr);
    ierr = VecSetFromOptions(Mnew);CHKERRQ(ierr);
    ierr = VecSetBlockSize(Mnew,ND);CHKERRQ(ierr);

    int cnt2;
    ierr = VecGetLocalSize(gdata->M,(PetscInt*)&cnt2);CHKERRQ(ierr);
    assert(cnt2==ND*gdata->n_vert);
    ierr = VecGetSize(gdata->M,(PetscInt*)&cnt2);CHKERRQ(ierr);
    assert(cnt2==ND*gdata->n_vert);

    ierr = VecGetLocalSize(Mnew,(PetscInt*)&cnt2);CHKERRQ(ierr);
    assert(cnt2==ND*(gdata->n_vert+n_edge));

    PetscReal *t_M, *t_Mnew;
    ierr = VecGetArray(gdata->M,&t_M);CHKERRQ(ierr);
    ierr = VecGetArray(Mnew,&t_Mnew);CHKERRQ(ierr);

    PetscPrintf(PETSC_COMM_WORLD,"Updating M\n");

    for (int i=0;i<gdata->n_vert;i++) {
      for (int k=0;k<ND;k++) {
        t_Mnew[ND*i+k]=t_M[ND*i+k];
      }
      for (int j=nxadj[i];j<nxadj[i+1];j++) {
        if (i<nadjncy[j]) {
          for (int k=0;k<ND;k++) {
            t_Mnew[ND*nvert[j]+k]=
              (t_M[ND*i+k]+t_M[ND*nadjncy[j]+k])/2.0;
          }
        }
      }
    }

    ierr = VecRestoreArray(gdata->M,&t_M);CHKERRQ(ierr);
    ierr = VecRestoreArray(Mnew,&t_Mnew);CHKERRQ(ierr);
    ierr = VecDestroy(gdata->M);CHKERRQ(ierr);
    gdata->M=Mnew;

    ierr = PetscFree(nvert);CHKERRQ(ierr);

    PetscReal devNorm;
    devNorm=RenormVec(gdata->M,1.0,0.0,ND);


    int *neleprop;
    ierr = PetscMalloc(8*gdata->n_ele*sizeof(int),&neleprop);CHKERRQ(ierr);
    for (int i=0;i<gdata->n_ele;i++) {
      for (int j=0;j<8;j++) {
        neleprop[8*i+j]=gdata->eleprop[i];
      }
    }
    ierr = PetscFree(gdata->eleprop);CHKERRQ(ierr);
    gdata->eleprop=neleprop;

    ierr = PetscFree(nxadj);CHKERRQ(ierr);
    ierr = PetscFree(nadjncy);CHKERRQ(ierr);

    /* update counters */
    gdata->n_vert  += n_edge;
    gdata->ln_vert = gdata->n_vert;

    gdata->n_ele   *= 8;
    gdata->ln_ele  = gdata->n_ele;
  }


  MagparFunctionLogReturn(0);
}