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