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

#include "field.h"

static int       doinit=1;
static int       fieldon=0;
static Vec       VHthis;
static PetscReal hextamp=0.5; /* Field strength in Tesla */

int Hext_cu_Init(GridData *gdata)
{
  PetscTruth     flg;

  MagparFunctionLogBegin;

  ierr = PetscOptionsGetInt(PETSC_NULL,"-hext_cu",(PetscInt*)&fieldon,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    fieldon=0;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -hext_cu not found, using default value: %i (field off)\n",
      fieldon
    );
  }
  if (!fieldon) {
    PetscPrintf(PETSC_COMM_WORLD,"Hext_cu off\n");
    MagparFunctionLogReturn(0);
  }
  PetscPrintf(PETSC_COMM_WORLD,"Hext_cu on\n");

  ierr = VecDuplicate(gdata->M,&VHthis);CHKERRQ(ierr);

  /* this example sets the external field in a certain area to a
     non-zero value and everywhere else to zero
  */

  /* make the vector of the external field accessible */
  PetscReal *ta_Hext;
  ierr = VecGetArray(VHthis,&ta_Hext);CHKERRQ(ierr);

  /* loop over all (local) nodes */
  for (int j=0; j<gdata->ln_vert; j++) {
    /* select nodes within a rectangular area */
    if (
      gdata->vertxyz[ND*gdata->vertl2g[j]+0]>=65.0 &&
      gdata->vertxyz[ND*gdata->vertl2g[j]+0]<=250.0 &&
      gdata->vertxyz[ND*gdata->vertl2g[j]+1]>=192.0 &&
      gdata->vertxyz[ND*gdata->vertl2g[j]+1]<=310.0
    ) {
      /* set the direction of the field (unit vector!) */

      PetscReal hextxyz[ND]={1.0,0.0,0.0};
      for (int i=0; i<ND; i++) {
        ta_Hext[ND*j+i]=hextxyz[i];
      }
    }
    else {
      /* set the field to zero outside the rectangular area */
      for (int i=0; i<ND; i++) {
        ta_Hext[ND*j+i]=0.0;
      }
    }
  }
  ierr = VecRestoreArray(VHthis,&ta_Hext);CHKERRQ(ierr);
  
  /* scale the vector to set the field strength */
  ierr = VecScale(gdata->VHext,hextamp);CHKERRQ(ierr);

  MagparFunctionLogReturn(0);
}


int Hext_cu(GridData *gdata,Vec VHextsum,PetscReal *hext)
{
  MagparFunctionInfoBegin;

  /* to be implemented by the user */

  if (doinit) {
    ierr = Hext_cu_Init(gdata);CHKERRQ(ierr);
    doinit=0;
  }
  if (!fieldon) {
    MagparFunctionInfoReturn(0);
  }

  ierr = VecAXPY(VHextsum,1.0,VHthis);CHKERRQ(ierr);

  *hext=hextamp;

  MagparFunctionProfReturn(0);
}