2nd order method

[Lin Ji]

Hi, Christophe and Jean-Francois,
    Sorry to bother you again.
    I tried the 2nd order interpolation and got some weird result. The
solution value becomes extremingly large (on the order of 10^7) around the
source as the wave propogates out. Can you take a look at the two problem
description files : 'Wave_u.pro' and 'test.pro' to see if I get it right? For
the constraint 'u2', I just copied the constraint 'u'.

Thanks,

Lin

-- 
Dept. of Math., RPI
518-276,2184(work)
518-271,4486(home)
-------------- next part --------------

/*
  A simple scalar wave propagation example
*/

Jacobian {
  { Name JVol ;
    Case { 
      { Region All ; Jacobian Vol ; }
    }
  }
  { Name JSur ;
    Case {
      { Region All ; Jacobian Sur ; }
    }
  }
}

Integration {
  { Name I1 ;
    Case { 
      { Type Gauss ;
        Case { 
	  { GeoElement Point       ; NumberOfPoints  1 ; }
	  { GeoElement Line        ; NumberOfPoints  3 ; }
	  { GeoElement Triangle    ; NumberOfPoints  4 ; }
	  { GeoElement Quadrangle  ; NumberOfPoints  4 ; }
	  { GeoElement Tetrahedron ; NumberOfPoints  4 ; }
	  { GeoElement Hexahedron  ; NumberOfPoints  6 ; }
	  { GeoElement Prism       ; NumberOfPoints  6 ; } 
	}
      }
    }
  }
}

FunctionSpace {
  { Name Hgrad; Type Form0;
    BasisFunction {
      { Name sn; NameOfCoef un; Function BF_Node; 
        Support Region[{Omega,Gamma_f}]; Entity NodesOf[All]; }
      { Name sn2; NameOfCoef un2; Function BF_Node_2E; 
        Support Region[{Omega,Gamma_f}]; Entity EdgesOf[All]; }
    }
    Constraint {
      { NameOfCoef un; EntityType NodesOf ; NameOfConstraint u; }
      { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint u2; }
   }
  }
}

Formulation {
  { Name Wave; Type FemEquation; 
    Quantity { 
      { Name u; Type Local; NameOfSpace Hgrad; }
    }
    Equation {
      Galerkin { DtDt [ 1/c[] * Dof{u} , {u} ]; 
                 In Omega; Integration I1; Jacobian JVol;  }

      Galerkin { [ c[] * Dof{d u} , {d u} ]; 
                 In Omega; Integration I1; Jacobian JVol;  }

      Galerkin {  [ - dfdt[] , {u} ]; 
                 In Gamma_f; Integration I1; Jacobian JSur;  }
    }
  }
}

Resolution {
  { Name Wave_t;
    System {
      { Name A; NameOfFormulation Wave; }
    }
    Operation { 
      InitSolution[A] ;
      InitSolution[A] ;
      TimeLoopNewmark[t_min,t_max,dt,0.25,0.5] {
	Generate[A] ; 
	Solve[A] ; 
	Test[ SaveFct[] ] {
	  SaveSolution[A] ; 
	}
      } 
    } 
  }


  { Name Wave_t_ex;
    System {
      { Name A; NameOfFormulation Wave; }
    }
    Operation { 
      InitSolution[A] ; 
      InitSolution[A] ;
      GenerateSeparate[A] ;
      TimeLoopNewmark[t_min,t_max,dt,0.25,0.5] {
	Update[A, TimeFct[]] ; 
	Solve[A] ; 
	If[ SaveFct[] ] {
	  SaveSolution[A] ;
	}
      } 
    } 
  }

  { Name Wave_f;
    System {
      { Name A; NameOfFormulation Wave; Type ComplexValue; Frequency Freq; }
    }
    Operation { 
      Generate[A]; 
      Solve[A];
      SaveSolution[A]; 
    }
  }

  { Name Wave_tf_ex;
    System {
      { Name A; NameOfFormulation Wave; }
      { Name C; NameOfFormulation Wave; Type ComplexValue ; }
    }
    Operation { 
      InitSolution[A] ;
      InitSolution[A] ;
      GenerateSeparate[A] ;
      TimeLoopNewmark[t_min,t_max,dt,0.25,0.5] {
	Update[A,TimeFct[]] ; 
	Solve[A] ; 
	FourierTransform[A, C, {Freq} ] ;
	If[SaveFct[]] { 
	  SaveSolution[A] ; 
	}
      } 
      SaveSolutions[C] ; 
    }    
  }

  { Name Wave_pvp;
    System {
      { Name A; NameOfFormulation Wave; }
    }
    Operation { 
      GenerateSeparate[A]; 
      Lanczos[A, 19, {1:5}, 0]; 
      SaveSolutions[A] ;
    }
  }
}


PostProcessing {
  { Name Wave_t; NameOfFormulation Wave; NameOfSystem A;
    Quantity {
      { Name u;  
        Value{ 
	  Local{ [ {u} ] ; In Omega; Jacobian JVol; } 
	} 
      }
      { Name du;  
        Value{ 
	  Local{ [ {d u} ] ; In Omega; Jacobian JVol; }
	} 
      }
      { Name u2; 
        Value{ 
	  Local{ [ Norm[{u}] ] ; In Omega; Jacobian JVol; } 
	} 
      }
    }
  }

  { Name Wave_f; NameOfFormulation Wave; NameOfSystem C;
    Quantity {
      { Name u;  
        Value{ 
	  Local{ [ {u} ] ; In Omega; Jacobian JVol; } 
	} 
      }
    }
  }
}


-------------- next part --------------

Group {
  Omega = Region[1] ;
  Gamma_u = Region[3] ;
  Gamma_f = Region[2] ;
}

Function {
 /*
 A=.5; 
 a=.01;
 c[] = 1+A*Exp[-((X[]-.5)^2+Y[]^2)/(2*a)] ;
 */
 c[] = 1;

  Freq = 1 ;  
  t_min = 0 ;
  t_max = .5 ;
  dt = .01 ;
  ct=.04;
  cy=.02;
  t0=.02;
  y0=0.0;
  gr=2*3.1416*100;
  yo[]=1/(ct*cy*Sqrt[2*3.1416])*Exp[-(($Time-t0)^2/(ct^2)+(Y[]-y0)^2/(cy^2) )];
  yo2[]=Cos[gr*($Time-t0)]*yo[];
  TimeFct[] = yo[];
  // TimeFct[] = ( (X[]<.03)  )? yo[] : 0 ;

  dfdt[] = TimeFct[] ;
 
  SaveFct[] = !($TimeStep % 2) ;
}

Constraint {
  { Name u ;
    Case {
      //{ Region Gamma_u  ; Value 0; }
      //{ Region Gamma_u  ; Value 1; TimeFunction TimeFct[]; }
    }
  }
  { Name u2 ;
    Case {
      //{ Region Gamma_u  ; Value 0; }
      //{ Region Gamma_u  ; Value 1; TimeFunction TimeFct[]; }
    }
  }
}

Include "Wave_u.pro"

PostOperation {
  { Name u ; NameOfPostProcessing Wave_t ;
    Operation {
      Print[ u, OnPlane { {0,-.5,0} {1,-.5,0} {0,.5,0} } {100,100} ,
       Format TimeTable,
       File "u.pos" ] ;
    }
  }

}

-------------- next part --------------
lc = 0.02;

Point(1) = {0,  -0.5,  0, lc};
Point(2) = {1, -0.5,  0, lc};
Point(3) = {1, .5, 0, lc};
Point(4) = {0,  .5, 0, lc};
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,1};
Line Loop(5) = {4,1,2,3};
Plane Surface(6) = {5};


Physical Surface(1) = {6};

Physical Line(2) = {4};
Physical Line(3) = {3};