Example ed4: 3D gradients from a potential domain solution

This example demonstrates potential domain solutions in three dimensions on a spherical BEM mesh solving Laplace's equation with an analytic solution. A perl subroutine is then used to calculate the potential gradients based on the potential field.

The comfile run by this example is as follows:

#
# Created by Martin Buist
# August 2001
#

use strict;

if( !$example ) {
  $example = ".";
}

fem define para;r;$example/3dtest
fem define coord 3,1          
fem define regi;r;$example/3dtest
fem define node;r;$example/8x8
fem define node;r;$example/3dtest reg 2
fem define node;w;3dtest reg 2
fem define bas;r;$example/test
fem define elem;r;$example/8x8
fem define elem;r;$example/3dtest reg 2
fem define grid;r;$example/3dtest reg 2
fem update grid geometry reg 2
fem define equation;r;$example/bemsphere cl 1 lock
fem define equa;r;$example/3dtest cl 2 lock reg 2
fem define material;r;$example/bemsphere 
fem define anal;r;$example/bemsphere
fem define initial;g 
fem define solve;r;$example/bemsphere
fem solve

fem eval solu grid elem 101 opti iy 1 potential

fem export node;geom as test reg 2
fem export elem;geom as test reg 2
fem export elem;pote as test field reg 2 class 2 iy 1

# Perl routines

sub ReadPotentialExelem {
 
  ( $INPUTNAME, $index, $xi1, $xi2, $xi3, $Nodes, $Values, @Potentials ) = @_;

  CORE::open( INPUTFILE, "$INPUTNAME" ) or die( "ERROR: Unable to open input file" );

  while( <INPUTFILE> ) {
    if( /xi1/ && /xi2/ ) {
      my @line = split;
      $xi1 = substr( $line[1], 0, -1 ) + 1;
      $xi2 = substr( $line[3], 0, -1 ) + 1;
      $xi3 = $line[5] + 1;
    }

    if( /Nodes:/ ) {
  		$Nodes = 0;
	  }

    if( $Values && $Nodes ) {
	  	my @line = split;
      foreach my $locindex ( 1..@line ) {  
	  		$Potentials[$index] = $line[$locindex-1];
  			$index++;
		  }
    } 

    if( /Values:/ ) {
  		$Values = 1;
	  }
  }

  CORE::close( INPUTFILE );

} # ReadPotentialExelem

sub CalcSteps {

  ( $NODENAME, $Xj1, $Xj2, $Xj3 ) = @_;
  my @XjExtremes = ( -10000, 10000, -10000, 10000, -10000, 10000 );

  CORE::open( INPUTFILE, "$NODENAME" ) or die( "ERROR: Unable to open input file" );

  while( <INPUTFILE> ) {
    # Minimum and maximum X coordinates
    if( /Xj\(1\)/ ) {
      my @line = split;
      $Xj1 = $line[-1];
      if ( $Xj1 > $XjExtremes[0] ) {
        $XjExtremes[0] = $Xj1;
      }
      if ( $Xj1 < $XjExtremes[1] ) {
        $XjExtremes[1] = $Xj1;
      }
    }

    # Minimum and maximum Y coordinates
    if( /Xj\(2\)/ ) {
      my @line = split;
      $Xj2 = $line[-1];
      if ( $Xj2 > $XjExtremes[2] ) {
        $XjExtremes[2] = $Xj2;
      }
      if ( $Xj2 < $XjExtremes[3] ) {
        $XjExtremes[3] = $Xj2;
      }
    }

    # Minimum and maximum Z coordinates
    if( /Xj\(3\)/ ) {
      my @line = split;
      $Xj3 = $line[-1];
      if ( $Xj3 > $XjExtremes[4] ) {
        $XjExtremes[4] = $Xj3;
      }
      if ( $Xj3 < $XjExtremes[5] ) {
        $XjExtremes[5] = $Xj3;
      }
    }
  }

  CORE::close( INPUTFILE );

  $Xj1 = ($XjExtremes[0] - $XjExtremes[1] / ($xi1 - 1));
  $Xj2 = ($XjExtremes[2] - $XjExtremes[3] / ($xi2 - 1));
  $Xj3 = ($XjExtremes[4] - $XjExtremes[5] / ($xi3 - 1));

} # CalcSteps

sub CalcGradPhi {

  ( $xi1, $xi2, $xi3, $POTE, $DX, $DY, $DZ ) = @_;
  my $Zero = 0.000001;

  foreach my $nik ( 1..$xi3 ) {
	  foreach my $nij ( 1..$xi2) {
		  foreach my $nii ( 1..$xi1 ) {
        $nq = $nii + (($nij-1)*$xi1) + (($nik-1)*$xi1*$xi2);
        $nqp  = $POTE->[$nq - 1];

        # Calculate dPhi/dx
        if( $nii == 1 ) {
          $nqp1 = $POTE->[$nq + 1 - 1];
          $nqp2 = $POTE->[$nq + 2 - 1];
          $DX->[$nq-1] = ((-3 * $nqp) + (4 * $nqp1) - $nqp2) / (2 * $Xj1);
        } elsif( $nii == $xi1 ) {
          $nqp1 = $POTE->[$nq - 1 - 1];
          $nqp2 = $POTE->[$nq - 2 - 1];
          $DX->[$nq-1] = ((3 * $nqp) - (4 * $nqp1) + $nqp2) / (2 * $Xj1);
        } else {
          $nqp1 = $POTE->[$nq + 1 - 1];
          $nqp2 = $POTE->[$nq - 1 - 1];
          $DX->[$nq-1] = ($nqp1 - $nqp2) / (2 * $Xj1);     
        }
        if( abs($nqp) <= $Zero || abs($nqp1) <= $Zero || abs($nqp2) <= $Zero ) {
          $DX->[$nq-1] = 0.0;
        }

        # Calculate dPhi/dy
        if( $nij == 1 ) {
          $nqp1 = $POTE->[$nq + $xi1 - 1];
          $nqp2 = $POTE->[$nq + (2 * $xi1) - 1];
          $DY->[$nq-1] = ((-3 * $nqp) + (4 * $nqp1) - $nqp2) / (2 * $Xj2);
        } elsif( $nij == $xi2 ) {
          $nqp1 = $POTE->[$nq - $xi1 - 1];
          $nqp2 = $POTE->[$nq - (2 * $xi1) - 1];
          $DY->[$nq-1] = ((3 * $nqp) - (4 * $nqp1) + $nqp2) / (2 * $Xj2);
        } else {
          $nqp1 = $POTE->[$nq + $xi1 - 1];
          $nqp2 = $POTE->[$nq - $xi1 - 1];
          $DY->[$nq-1] = ($nqp1 - $nqp2) / (2 * $Xj2);  
        }
        if( abs($nqp) <= $Zero || abs($nqp1) <= $Zero || abs($nqp2) <= $Zero ) {
          $DY->[$nq-1] = 0.0;
        }

        # Calculate dPhi/dz
        if( $nik == 1 ) {
          $nqp1 = $POTE->[$nq + ($xi1 * $xi2) - 1];
          $nqp2 = $POTE->[$nq + (2 * $xi1 * $xi2) - 1];
          $DZ->[$nq-1] = ((-3 * $nqp) + (4 * $nqp1) - $nqp2) / (2 * $Xj3);
        } elsif( $nik == $xi3 ) {
          $nqp1 = $POTE->[$nq + ($xi1 * $xi2) - 1];
          $nqp2 = $POTE->[$nq - (2 * $xi1 * $xi2) - 1];
          $DZ->[$nq-1] = ((3 * $nqp) - (4 * $nqp1) + $nqp2) / (2 * $Xj3);
        } else {
          $nqp1 = $POTE->[$nq + ($xi1 * $xi2) - 1];
          $nqp2 = $POTE->[$nq - ($xi1 * $xi2) - 1];
          $DZ->[$nq-1] = ($nqp1 - $nqp2) / (2 * $Xj3);    
        }
        if( abs($nqp) <= $Zero || abs($nqp1) <= $Zero || abs($nqp2) <= $Zero ) {
          $DZ->[$nq-1] = 0.0;
        }
      }
    }
  }

} # CalcGradPhi

sub WriteGradientExelem {

  ( $Fields, $Values, $INPUTNAME, $OUTPUTNAME, $xi1, $xi2, $xi3, $index, $DX, $DY, $DZ ) = @_;
  $xi1--;
  $xi2--;
  $xi3--;
  $Values = 0;
 
  CORE::open( INPUTFILE, "$INPUTNAME" ) or die( "ERROR: Unable to open input file" );
  CORE::open( OUTPUTFILE, ">$OUTPUTNAME" ) or die( "ERROR: Unable to open output file" );

  while( <INPUTFILE> ) {
    if( /1\)/ ) {
      $Fields = 1;
      print( OUTPUTFILE ' 1) gradphi, field, rectangular cartesian, #Components=3'."\n" );
      print( OUTPUTFILE '   dx.  l.Lagrange*l.Lagrange*l.Lagrange, no modify, grid based.'."\n" );
      print( OUTPUTFILE '   #xi1=  '."$xi1".', #xi2=  '."$xi2".', #xi3= '."$xi3 \n" );
      print( OUTPUTFILE '   dy.  l.Lagrange*l.Lagrange*l.Lagrange, no modify, grid based.'."\n" );
      print( OUTPUTFILE '   #xi1=  '."$xi1".', #xi2=  '."$xi2".', #xi3= '."$xi3 \n" );
      print( OUTPUTFILE '   dz.  l.Lagrange*l.Lagrange*l.Lagrange, no modify, grid based.'."\n" );
      print( OUTPUTFILE '   #xi1=  '."$xi1".', #xi2=  '."$xi2".', #xi3= '."$xi3 \n" );
  	}

    if( /Element:/ ) {
      $Fields = 0;
    }

    if( /Nodes:/ ) {
      $Values = 0;
    }
  
    if( /Values:/ ) {
      $Values = 1;
      print( OUTPUTFILE "   Values:\n    " );
    
      foreach my $locindex ( 1..$index ) {
        print( OUTPUTFILE "  $DX->[$locindex-1]" );
        if( $locindex % 5 == 0 ) {
          print( OUTPUTFILE "\n    " );
        }
      }
      foreach my $locindex ( 1..$index ) {
        print( OUTPUTFILE "  $DY->[$locindex-1]" );
        if( $locindex % 5 == 0 ) {
          print( OUTPUTFILE "\n    " );
        }
      }
      foreach my $locindex ( 1..$index ) {
        print( OUTPUTFILE "  $DZ->[$locindex-1]" );
        if( $locindex % 5 == 0 ) {
          print( OUTPUTFILE "\n    " );
        }
      }  
    }

    if( !$Fields && !$Values ) {
      print( OUTPUTFILE $_  );
    }
  }

  CORE::close( INPUTFILE );
  CORE::close( OUTPUTFILE );

} # WriteGradientExelem

# Main Loop

$index = 0;
$xi1 = 0;
$xi2 = 0;
$xi3 = 0;
$Xj1 = 0;
$Xj2 = 0;
$Xj3 = 0;
@Potentials = ' ';
$Values = 0;
$Nodes = 1;
$INPUTNAME = "pote.exelem";
$NODENAME = "3dtest.ipnode";
$OUTPUTNAME = "FDgradphi.exelem";
$Fields = 0;
@dPhidX = ' ';
@dPhidY = ' ';
@dPhidZ = ' ';

&ReadPotentialExelem( $INPUTNAME, $index, $xi1, $xi2, $xi3, $Nodes, $Values, @Potentials );
&CalcSteps( $NODENAME, $Xj1, $Xj2, $Xj3 );
$DX = \@dPhidX;
$DY = \@dPhidY;
$DZ = \@dPhidZ;
$POTE = \@Potentials
&CalcGradPhi( $xi1, $xi2, $xi3, $POTE, $DX, $DY, $DZ );
&WriteGradientExelem( $Fields, $Values, $INPUTNAME, $OUTPUTNAME, $xi1, $xi2, $xi3, $index, $DX, $DY, $DZ );

Additional testing commands:








Files used by this example are:





Name              Modified     Size


example_ed4.com   04-Sep-2001  7.8k
3dtest.ipbase     04-Sep-2001  2.4k
3dtest.ipelem     04-Sep-2001  422 
3dtest.ipequa     26-May-2003  1.2k
3dtest.ipgrid     06-Mar-2003  605 
3dtest.ipnode     04-Sep-2001  2.2k
3dtest.ippara     12-Nov-2002  5.9k
3dtest.ipregi     04-Sep-2001   93 
8x8.ipelem        04-Sep-2001   21k
8x8.ipnode        04-Sep-2001   52k
FDgradphi.exelem  04-Sep-2001   20k
bemsphere.ipanal  04-Sep-2001  591 
bemsphere.ipequa  26-May-2003  1.4k
bemsphere.ipmate  04-Sep-2001   54 
bemsphere.ipsolv  13-Apr-2007  1.3k
geom.exelem       04-Sep-2001  3.6k
geom.exnode       04-Sep-2001  1.0k
pote.exelem       04-Sep-2001   15k
test.ipbase       04-Sep-2001  6.8k
test_output.com   04-Sep-2001    0 








Download the entire example:





Name                      Modified     Size


examples_e_ed_ed4.tar.gz  16-Aug-2014   19k









Testing status by version:




 StatusTestedReal time (s)


i686-linux


cmSuccessSun Mar  6 00:04:43 20165


cm-debugSuccessSat Mar  5 00:25:17 201618


mips-irix


cmSuccessSun Aug 19 02:13:49 200743


cm-debugSuccessWed Aug 15 02:46:29 2007195


cm-debug-clear-mallocSuccessSat Aug 18 03:02:20 2007197


cm-debug-clear-malloc7SuccessMon Aug 20 03:03:10 2007212


cm64SuccessSun Aug 19 02:19:29 200747


cm64-debugSuccessTue Aug 21 03:10:17 2007215


cm64-debug-clear-mallocSuccessThu Apr  1 12:31:24 2004118


rs6000-aix


cmSuccessWed Mar  4 01:12:13 20097


cm-debugSuccessMon Mar  2 01:24:12 200958


cm64SuccessWed Mar  4 01:12:32 20096


cm64-debugSuccessTue Mar  3 01:27:37 200962


x86_64-linux


cmSuccessSun Mar  6 00:01:13 20162


cm-debugSuccessSat Mar  5 00:02:20 20169






Testing status by file:








 Html last generated: Sun Mar  6 05:51:36 2016



 Input last modified: Fri Aug 15 13:13:24 2014






CMISS Help /
Examples /
e /
ed /
ed4