Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: how to integrate a vector over a discrete surface?
Date: Tue, 8 Mar 2011 21:48:19 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 36
Message-ID: <il6873$mhm$1@fred.mathworks.com>
References: <il3k7n$dqu$1@fred.mathworks.com> <il459e$opn$1@fred.mathworks.com> <il5gld$2oh$1@fred.mathworks.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: www-02-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1299620899 23094 172.30.248.47 (8 Mar 2011 21:48:19 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Tue, 8 Mar 2011 21:48:19 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: news.mathworks.com comp.soft-sys.matlab:714663

"yl " <shanshicn@hotmail.com> wrote in message <il5gld$2oh$1@fred.mathworks.com>...
> Roger,Thanks a lot. That is very helpful. Just to understand better
> about the principle, any recommend readings? thanks.

  I'm not sure what reading to recommend to you.  (My own sources are very ancient.)  I'll show you an intuitive derivation of the formula I gave and that should at least be an indication of where in more modern literature to look.

  The line integral formula "F = P/2*int (r cross dr)" is a consequence of the famous Green's Theorem.  See:

 http://en.wikipedia.org/wiki/Green's_theorem

  If S is the surface you are integrating the constant pressure, P, over, let ds be an infinitesimal normal vector area element.  Then the orthogonal projection of P*ds onto the x-y plane would just be the infinitesimal z-component of the force exerted on that surface element.  Therefore the total z-component of the force would simply be the projected area of S onto the x-y plane times P which is bounded by the similar projection of the S boundary there.  (If the projected surface bulges past this projected boundary there will always be an appropriate cancellation of that portion, so we only need consider the interior of the projected boundary.)

  The z-component of the formula's "P/2*int (r cross dr)" will be

 P/2*int (x*dy-y*dx)

integrated along the projected boundary.  By Green's Theorem using their L = -y and M = x (see the above website) this would be equal to

 P/2*int (dM/dx-dL/dy) dxdy = P/2*int (1+1) dxdy = P*int 1 dxdy

integrated over the projected x-y area, which will give us just the projected (interior) area times P, and therefore the z-component of the total force.

  Applied to all three x, y, and z components, this justifies the above formula.

> Roger, one more question. which integration command in matlab for the
> integration on discrete boundary? thanks a lot.

  That question is not so easy to answer.  In your original x,y,z data, to do a proper surface integration you would have had to know accurately where the boundary of the surface lay.  With your data at scattered interior discrete points, it would not be evident just where this was.

  Presuming you somehow obtain a sequence of sufficiently closely-spaced successive points lying right along this required boundary, you could use the 'trapz' function for each of the three components of P/2*int (r cross dr) to find the three force components.  The remaining data points of the interior would not be needed at all.  If such boundary points were not sufficiently close together to achieve the desired accuracy, you might have to resort to some interpolation to provide more closely-spaced points for 'trapz' or perhaps you could use one of the higher order integration schemes for discrete points present on the file exchange.  I wrote one at:

 http://www.mathworks.com/matlabcentral/fileexchange/19152

which does third-order cumulative integration of a sequence of discrete points, and there are a number of others to be found there.

Roger Stafford