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

#include "init.h"
#include "io/magpario.h"

int ModifyPropPar(GridData *gdata)
{
  MagparFunctionLogBegin;

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

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

  /* exit if there is nothing to modify */
  if (nsliceprop<=0) {
    PetscPrintf(PETSC_COMM_WORLD,
      "nsliceproppar==%i, nothing modified!\n",
      nsliceprop
    );
    MagparFunctionLogReturn(0);
  }

  /* sliced4.out */
  PetscReal ep0, ew; /* element "position", width */
  ep0= 0.1;
  ew=  0.25;

  int ps0, ns; /* prop. id of first slice, no. of slices (properties) */
  ps0=15;
  ns= 40;

  PetscPrintf(PETSC_COMM_WORLD,
    "modifyprop: element position: %g width: %g propid: [%i,%i]\n",
    ep0,
    ew,
    ps0,
    ps0+ns-1
  );

  if (ps0<gdata->n_prop) {
    SETERRQ1(PETSC_ERR_ARG_INCOMP,"ps0 < gdata->n_prop in %s\n",__FUNCT__);
  }

  if (ns != nsliceprop) {
    SETERRQ1(PETSC_ERR_ARG_INCOMP,"ns != nsliceproppar in %s\n",__FUNCT__);
  }

  gdata->n_prop += ns;

  PetscPrintf(PETSC_COMM_WORLD,
    "modifyprop: rereading %s.krn for %i additional properties\n",
    gdata->simname,
    ns
  );
  ierr = PetscFree(gdata->propdat);CHKERRQ(ierr);
  ierr = ReadKrn(gdata);CHKERRQ(ierr);

  if (rank) {
    ierr = PetscMalloc(gdata->n_prop*NP*sizeof(PetscReal),&gdata->propdat);CHKERRQ(ierr);
  }
  ierr = MPI_Bcast(gdata->propdat,gdata->n_prop*NP,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);

  for (int j=0; j<ns; j++) {
    for (int i=0; i<gdata->ln_ele; i++) {
      int k;
      /* get x-coordinate of first vertex */
      PetscReal sgn;
      sgn=gdata->vertxyz[ND*gdata->elevert[NV*i+0]+0]-(ep0+ew*j);
      for (k=1;k<NV;k++) {
        /* multiply with other x-coordinate of other vertices
           if (product < 0) then they have opposite signs and
           the element must be cut by the slice plane
        */
        if (sgn*(gdata->vertxyz[ND*gdata->elevert[NV*i+k]+0]-(ep0+ew*j)) < 0.0)
          break;
      }
      if (k<NV) {
        gdata->eleprop[i]=ps0+j;
      }

/* does not always work (due to strongly deformed elements)
      if (PetscAbsReal(
            PetscAbsReal(gdata->vertxyz[ND*gdata->elevert[NV*i+0]+0]+
                 gdata->vertxyz[ND*gdata->elevert[NV*i+1]+0]+
                 gdata->vertxyz[ND*gdata->elevert[NV*i+2]+0]+
                 gdata->vertxyz[ND*gdata->elevert[NV*i+3]+0]
                )/NV-(ep-j*ew)
              ) <= et
         ) {
        gdata->eleprop[i]=ps0+j;
      }
*/
    }
  }

  MagparFunctionLogReturn(0);
}