How to code a nonlinear vector-valued function so that it fits the coordinate grids (xq , yq) of triangles?

This code got me chalenged. My problem is how to re-code the function fq = f ( xq, yq ) into a nonlinear function fq = f ( u (xq, yq) ) so that it fits the coordinate grids (xq , yq) of triangular elements? Also, If the boundary condition is a nonzero function, what to change in the code?

function u = fem2d_bvp_linear ( nx, ny, a, c, f, x, y )

%*****************************************************************************80

%

%% fem2d_bvp_linear() solves boundary value problem on a rectangle.

%

% Discussion:

%

% The program uses the finite element method, with piecewise linear basis

% functions to solve a 2D boundary value problem over a rectangle

%

% The following differential equation is imposed inside the region:

%

% - d/dx a(x,y) du/dx - d/dy a(x,y) du/dy + c(x,y) * u(x,y) = f(x,y)

%

% where a(x,y), c(x,y), and f(x,y) are given functions.

%

% On the boundary, the solution is constrained to have the value 0.

%

% The finite element method will use a regular grid of NX nodes in X, and

% NY nodes in Y.

%

% Licensing:

%

% This code is distributed under the GNU LGPL license.

%

% Modified:

%

% 20 June 2014

%

% Author:

%

% John Burkardt

%

% Input:

%

% integer NX, NY, the number of X and Y grid values.

%

% function A(X,Y), evaluates a(x,y);

%

% function C(X,Y), evaluates c(x,y);

%

% function F(X,Y), evaluates f(x,y);

%

% real X(NX), Y(NY), the mesh points.

%

% Output:

%

% real U(NX,NY), the finite element coefficients, which are also

% the value of the computed solution at the mesh points.

%

%

% Quadrature definitions.

%

quad_num = 3;

abscissa(1) = -0.774596669241483377035853079956;

abscissa(2) = 0.000000000000000000000000000000;

abscissa(3) = 0.774596669241483377035853079956;

weight(1) = 0.555555555555555555555555555556;

weight(2) = 0.888888888888888888888888888889;

weight(3) = 0.555555555555555555555555555556;

%

% Make room for the matrix A and right hand side b.

%

mn = nx * ny;

A = zeros ( mn, mn );

b = zeros ( mn, 1 );

%

% Compute the matrix entries by integrating over each element.

%

for ex = 1 : nx - 1

w = ex;

e = ex + 1;

xw = x(w);

xe = x(e);

for ey = 1 : ny - 1

s = ey;

n = ey + 1;

ys = y(s);

yn = y(n);

sw = ( ey - 1 ) * nx + ex;

se = ( ey - 1 ) * nx + ex + 1;

nw = ey * nx + ex;

ne = ey * nx + ex + 1;

for qx = 1 : quad_num

xq = ( ( 1.0 - abscissa(qx) ) * xw ...

+ ( 1.0 + abscissa(qx) ) * xe ) ...

/ 2.0;

for qy = 1 : quad_num

yq = ( ( 1.0 - abscissa(qy) ) * ys ...

+ ( 1.0 + abscissa(qy) ) * yn ) ...

/ 2.0;

wq = weight(qx) * ( xe - xw ) / 2.0 ...

* weight(qy) * ( yn - ys ) / 2.0;

vsw = ( xe - xq ) / ( xe - xw ) * ( yn - yq ) / ( yn - ys );

vswx = ( -1.0 ) / ( xe - xw ) * ( yn - yq ) / ( yn - ys );

vswy = ( xe - xq ) / ( xe - xw ) * ( -1.0 ) / ( yn - ys );

vse = ( xq - xw ) / ( xe - xw ) * ( yn - yq ) / ( yn - ys );

vsex = ( 1.0 ) / ( xe - xw ) * ( yn - yq ) / ( yn - ys );

vsey = ( xq - xw ) / ( xe - xw ) * ( -1.0 ) / ( yn - ys );

vnw = ( xe - xq ) / ( xe - xw ) * ( yq - ys ) / ( yn - ys );

vnwx = ( -1.0 ) / ( xe - xw ) * ( yq - ys ) / ( yn - ys );

vnwy = ( xe - xq ) / ( xe - xw ) * ( 1.0 ) / ( yn - ys );

vne = ( xq - xw ) / ( xe - xw ) * ( yq - ys ) / ( yn - ys );

vnex = ( 1.0 ) / ( xe - xw ) * ( yq - ys ) / ( yn - ys );

vney = ( xq - xw ) / ( xe - xw ) * ( 1.0 ) / ( yn - ys );

aq = a ( xq, yq );

cq = c ( xq, yq );

% fq = f ( xq, yq ); How to code this if I have fq = f(u(xq,yq))

% where u is the solution?

A(sw,sw) = A(sw,sw) + wq * ( vswx * aq * vswx ...

+ vswy * aq * vswy ...

+ vsw * cq * vsw );

A(sw,se) = A(sw,se) + wq * ( vswx * aq * vsex ...

+ vswy * aq * vsey ...

+ vsw * cq * vse );

A(sw,nw) = A(sw,nw) + wq * ( vswx * aq * vnwx ...

+ vswy * aq * vnwy ...

+ vsw * cq * vnw );

A(sw,ne) = A(sw,ne) + wq * ( vswx * aq * vnex ...

+ vswy * aq * vney ...

+ vsw * cq * vne );

b(sw) = b(sw) + wq * ( vsw * fq );

A(se,sw) = A(se,sw) + wq * ( vsex * aq * vswx ...

+ vsey * aq * vswy ...

+ vse * cq * vsw );

A(se,se) = A(se,se) + wq * ( vsex * aq * vsex ...

+ vsey * aq * vsey ...

+ vse * cq * vse );

A(se,nw) = A(se,nw) + wq * ( vsex * aq * vnwx ...

+ vsey * aq * vnwy ...

+ vse * cq * vnw );

A(se,ne) = A(se,ne) + wq * ( vsex * aq * vnex ...

+ vsey * aq * vney ...

+ vse * cq * vne );

b(se) = b(se) + wq * ( vse * fq );

A(nw,sw) = A(nw,sw) + wq * ( vnwx * aq * vswx ...

+ vnwy * aq * vswy ...

+ vnw * cq * vsw );

A(nw,se) = A(nw,se) + wq * ( vnwx * aq * vsex ...

+ vnwy * aq * vsey ...

+ vnw * cq * vse );

A(nw,nw) = A(nw,nw) + wq * ( vnwx * aq * vnwx ...

+ vnwy * aq * vnwy ...

+ vnw * cq * vnw );

A(nw,ne) = A(nw,ne) + wq * ( vnwx * aq * vnex ...

+ vnwy * aq * vney ...

+ vnw * cq * vne );

b(nw) = b(nw) + wq * ( vnw * fq );

A(ne,sw) = A(ne,sw) + wq * ( vnex * aq * vswx ...

+ vney * aq * vswy ...

+ vne * cq * vsw );

A(ne,se) = A(ne,se) + wq * ( vnex * aq * vsex ...

+ vney * aq * vsey ...

+ vne * cq * vse );

A(ne,nw) = A(ne,nw) + wq * ( vnex * aq * vnwx ...

+ vney * aq * vnwy ...

+ vne * cq * vnw );

A(ne,ne) = A(ne,ne) + wq * ( vnex * aq * vnex ...

+ vney * aq * vney ...

+ vne * cq * vne );

b(ne) = b(ne) + wq * ( vne * fq );

end

end

end

end

%

% Where a node is on the boundary,

% replace the finite element equation by a boundary condition.

%

k = 0;

for y = 1 : ny

for x = 1 : nx

k = k + 1;

if ( x == 1 || x == nx || y == 1 || y == ny )

A(k,1:mn) = 0.0;

A(1:mn,k) = 0.0;

A(k,k) = 1.0;

b(k) = 0.0;

end

end

end

if ( false )

spy ( A );

pause

end

%

% Solve the linear system.

%

u = A \ b;

Answer 1

[Torsten]

Here is the author's homepage:

https://people.sc.fsu.edu/~jburkardt/

Just contact him.

[Qais Al Faraei]

Fine.. Thanks a lot