l2-rhs-f-grad-g-cartesian-q2.at

这是FreeFem2D的姊妹软件

# This file is part of ff3d - http://www.freefem.org/ff3d
# Copyright (C) 2005 Stephane Del Pino

# This program 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, or (at your option)
# any later version.

# This program 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 this program; if not, write to the Free Software Foundation,
# Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

# $Id: l2-rhs-f-grad-g-cartesian-q2.at,v 1.1 2007/06/19 22:53:11 delpinux Exp $

AT_SETUP([L2 rhs: f*grad g          (cartesian Q2-FEM)])

AT_DATA([test.ff],
[[
vertex a = (1,0,0);
vertex b = (3,2,3);
vertex n = (3,2,3);

mesh M = structured(n,a,b);

double a1 = 1.5;

double a2 = 01;
double b2 = 02;
double c2 = 03;
double d2 = 04;
double e2 = 05;
double f2 = 06;
double g2 = 07;

double a3 = 11;
double b3 = 12;
double c3 = 13;
double d3 = 14;
double e3 = 15;
double f3 = 16;
double g3 = 17;

double a4 = 21;
double b4 = 22;
double c4 = 23;
double e4 = 25;

double a5 = 31;
double b5 = 32;
double c5 = 33;
double d5 = 34;
double e5 = 35;
double f5 = 36;
double g5 = 37;

double a6 = 41;
double b6 = 42;
double c6 = 43;
double f6 = 46;

function fx = a1*x;
function gx = a2*x + b2*y + c2*z + d2*y*z + e2*x*z + f2*x*y + g2*x*y*z;

function fy = a3*x + b3*y + c3*z + d3*y*z + e3*x*z + f3*x*y + g3*x*y*z;
function gy = a4*x + b4*y + c4*z +          e4*x*z;

function fz = a5*x + b5*y + c5*z + d5*y*z + e5*x*z + f5*x*y + g5*x*y*z;
function gz = a6*x + b6*y + c6*z + f6*x*y;

function uexact = (a2+e2*z+f2*y+g2*y*z)*fx + b4*fy + c6*fz;

solve(u:P2) in M
  cg(epsilon=1E-30)
{
  test(v)
   int(u*v) = int(fx*dx(gx)*v + fy*dy(gy)*v  + fz*dz(gz)*v);
   u =  uexact on M;
}

double l2error = sqrt(int[M:P2]((u-uexact)^2));
if (l2error < 1E-8)
  cout << "ok\n";
else
  cout << "oops: l2 error is " << l2error << "\n";
]])

AT_CHECK([$abs_top_builddir/ff3d -V 0 -nw test.ff],0,
[[ok
]])

AT_CLEANUP