distortmesh.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: distortmesh.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "init.h"

int DistortMesh(GridData *gdata)
{
  PetscTruth   flg;

  MagparFunctionLogBegin;

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

  if (rank) {
    MagparFunctionLogReturn(0);
  }

  /* scale mesh */

  PetscReal mesh_scale[ND];
  int nmax;

  nmax=ND;
  ierr = PetscOptionsGetRealArray(PETSC_NULL,"-mesh_scale",mesh_scale,&nmax,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    /* set default values if option is missing */
    mesh_scale[0]=mesh_scale[1]=mesh_scale[2]=1.0;
    PetscPrintf(PETSC_COMM_WORLD,"Option -mesh_scale not found or incorrect number of values given.\nUsing default value '%g,%g,%g'!\n",
      mesh_scale[0],mesh_scale[1],mesh_scale[2]
    );
  }
  else if (nmax!=ND) {
    /* set default values if option is missing */
    SETERRQ1(PETSC_ERR_ARG_CORRUPT,
      "Option -mesh_scale has incorrect number of values: %i!\n",
      nmax
    );
  }
  else {
    PetscPrintf(PETSC_COMM_WORLD,"Scaling mesh: '%g,%g,%g'!\n",
      mesh_scale[0],mesh_scale[1],mesh_scale[2]
    );
    for (int i=0; i<gdata->ln_vert; i++) {
      for (int j=0; j<ND; j++) {
        gdata->vertxyz[ND*i+j] *= mesh_scale[j];
      }
    }
  }


  /* shift mesh */

  PetscReal shift[ND];
  nmax=ND;
  ierr = PetscOptionsGetRealArray(PETSC_NULL,"-shift",shift,&nmax,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    /* set default values if option is missing */
    shift[0]=shift[1]=shift[2]=0;
    PetscPrintf(PETSC_COMM_WORLD,"Option -shift not found, using default value '%g,%g,%g'!\n",
      shift[0],shift[1],shift[2]
    );
  }
  else if (nmax!=ND) {
    /* set default values if option is missing */
    SETERRQ1(PETSC_ERR_ARG_CORRUPT,
      "Option -shift has incorrect number of values: %i!\n",
      nmax
    );
  }
  else {
    PetscPrintf(PETSC_COMM_WORLD,"Shifting mesh: '%g,%g,%g'!\n",
      shift[0],shift[1],shift[2]
    );
    for (int i=0; i<gdata->ln_vert; i++) {
      for (int j=0; j<ND; j++) {
        gdata->vertxyz[ND*i+j] += shift[j];
      }
    }
  }


  /* distort the mesh */

  int meshdist;
  ierr = PetscOptionsGetInt(PETSC_NULL,"-meshdist",(PetscInt*)&meshdist,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    /* set default values if option is missing */
    meshdist=0;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -meshdist not found, using default value %i!\n",
      meshdist
    );
  }

  PetscPrintf(PETSC_COMM_WORLD,"meshdist: %i\n",meshdist);
  if (meshdist==0) {
    PetscPrintf(PETSC_COMM_WORLD,"mesh undistorted\n");
    MagparFunctionLogReturn(0);
  }

  PetscReal distpar;
  ierr = PetscOptionsGetReal(PETSC_NULL,"-distpar",&distpar,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    /* set default values if option is missing */
    distpar=0;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -distpar not found, using default value %g!\n",
      distpar
    );
  }

  /* elevertvol will check if the mesh is corrupted
     i.e. overlapping elements = negative volumes of tetrahedra
  */

  PetscReal elenmax,elenmin;
  elenmax=0.0;
  elenmin=1.0/D_EPS;

  for (int i=0; i<gdata->ln_ele; i++) {
    /* Get Vertex coordinates */
    int *t_v;
    t_v=gdata->elevert+i*NV;

    PetscReal *x1,*x2,*x3,*x4;
    x1=gdata->vertxyz+t_v[0]*ND;
    x2=gdata->vertxyz+t_v[1]*ND;
    x3=gdata->vertxyz+t_v[2]*ND;
    x4=gdata->vertxyz+t_v[3]*ND;

    /* Calculate edge lengths */
    PetscReal edge[ND];
    my_dcopy(ND,x1,1,edge,1);
    my_daxpy(ND,-1.0,x2,1,edge,1);

    PetscReal elen;
    elen=my_dnrm2(ND,edge,1);
    elenmax=PetscMax(elenmax,elen);
    elenmin=PetscMin(elenmin,elen);

    my_dcopy(ND,x1,1,edge,1);
    my_daxpy(ND,-1.0,x3,1,edge,1);
    elen=my_dnrm2(ND,edge,1);
    elenmax=PetscMax(elenmax,elen);
    elenmin=PetscMin(elenmin,elen);

    my_dcopy(ND,x1,1,edge,1);
    my_daxpy(ND,-1.0,x4,1,edge,1);
    elen=my_dnrm2(ND,edge,1);
    elenmax=PetscMax(elenmax,elen);
    elenmin=PetscMin(elenmin,elen);

    my_dcopy(ND,x2,1,edge,1);
    my_daxpy(ND,-1.0,x3,1,edge,1);
    elen=my_dnrm2(ND,edge,1);
    elenmax=PetscMax(elenmax,elen);
    elenmin=PetscMin(elenmin,elen);

    my_dcopy(ND,x2,1,edge,1);
    my_daxpy(ND,-1.0,x4,1,edge,1);
    elen=my_dnrm2(ND,edge,1);
    elenmax=PetscMax(elenmax,elen);
    elenmin=PetscMin(elenmin,elen);

    my_dcopy(ND,x3,1,edge,1);
    my_daxpy(ND,-1.0,x4,1,edge,1);
    elen=my_dnrm2(ND,edge,1);
    elenmax=PetscMax(elenmax,elen);
    elenmin=PetscMin(elenmin,elen);
  }

  PetscPrintf(PETSC_COMM_WORLD,"distorting mesh - meshdist: %i distpar: %g elenmin: %g\n",
    meshdist,distpar,elenmin
  );

  PetscRandom rctx;
  ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
  ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);

  for (int i=0; i<gdata->ln_vert; i++) {
    PetscReal randist[ND];
    PetscReal norm;

    /* calculate random vector for M, which lies within
       the |M|=1 sphere, to guarantee uniform distribution
       on the unit sphere
    */
    do {
      PetscRandomGetValue(rctx,&randist[0]); randist[0] = randist[0]*2.0-1.0;
      PetscRandomGetValue(rctx,&randist[1]); randist[1] = randist[1]*2.0-1.0;
      PetscRandomGetValue(rctx,&randist[2]); randist[2] = randist[2]*2.0-1.0;
    } while ((norm=my_dnrm2(ND,randist,1)) > 1.0);

    switch (meshdist) {
      case 1: /* distort only interior mesh */
        if (gdata->vertbndg2bnd[gdata->vertl2g[i]] == C_INT) {
          gdata->vertxyz[ND*i+0] += randist[0]/norm*distpar*elenmin;
          gdata->vertxyz[ND*i+1] += randist[1]/norm*distpar*elenmin;
          gdata->vertxyz[ND*i+2] += randist[2]/norm*distpar*elenmin;
        }
      break;
      case 2: /* distort only boundary mesh */
        if (gdata->vertbndg2bnd[gdata->vertl2g[i]] >= 0) {
          gdata->vertxyz[ND*i+0] += randist[0]/norm*distpar*elenmin;
          gdata->vertxyz[ND*i+1] += randist[1]/norm*distpar*elenmin;
          gdata->vertxyz[ND*i+2] += randist[2]/norm*distpar*elenmin;
        }
      break;
      case 3: /* distort whole mesh */
        gdata->vertxyz[ND*i+0] += randist[0]/norm*distpar*elenmin;
        gdata->vertxyz[ND*i+1] += randist[1]/norm*distpar*elenmin;
        gdata->vertxyz[ND*i+2] += randist[2]/norm*distpar*elenmin;
      break;
    }
  }
  ierr = PetscRandomDestroy(rctx);CHKERRQ(ierr);

  MagparFunctionLogReturn(0);
}