| Author | Message | ||
     Matthew Hills (mdhills) Member Username: mdhills Post Number: 7 Registered: 01-2007 |
I'd like to calculate the surface area of a region that is defined in cylindrical coordinates. I thought that a simple surface integral could calculate this: report( SURF_INTEGRAL(1,"cathode")) as 'Cathode area (mm2)' However, when run, the result seems to depend on which direction I trace the boundary of the region -- is that right? In particular, if I trace the boundary clockwise, the area calcualtion ends up double that if it is done counter-clockwise. I also noticed that pulling the shapes off of the axis of symmetry by a small amount results in both integrals coming out pretty much identical... Am I missing something? E.g., { test_area.pde } TITLE 'Calculate areas' COORDINATES xcylinder(Z,R) SELECT alias(z) = "Z(mm)" alias(r) = "r(mm)" PAINTED VARIABLES { system variables } V DEFINITIONS rad = 0.5 EQUATIONS DIV(GRAD(V)) = 0 BOUNDARIES REGION 1 start ( -10, 0) natural(V) = 0 line to (10,0) value(V) = 0 line to (10,10) to (-10,10) to close REGION 2 start 'cw' (-3*rad,0) value(V)=-10 arc(center=-2*rad,0) to (-2*rad,rad) arc(center=-2*rad,0) to (-rad,0) natural(V)=0 line to close REGION 3 start 'ccw' (3*rad,0) value(V)=-10 arc(center=2*rad,0) to (2*rad,rad) arc(center=2*rad,0) to (rad,0) natural(V)=0 line to close PLOTS contour(V) zoom(-1,0, 2,3) SUMMARY report( SURF_INTEGRAL(1,"cw")) as 'Clockwise Area (mm2)' report( SURF_INTEGRAL(1,"ccw")) as 'Counter-clockwise Area (mm2)' END | ||
| Matthew Hills (mdhills) Member Username: mdhills Post Number: 8 Registered: 01-2007 |
Here's a similar example, showing a boundary that has been defined clockwise, but with a slight displacement from the axis of symmetry: { test_area.pde } TITLE 'Calculate areas' COORDINATES xcylinder(Z,R) SELECT alias(z) = "Z(mm)" alias(r) = "r(mm)" PAINTED VARIABLES { system variables } V DEFINITIONS rad = 0.5 EQUATIONS DIV(GRAD(V)) = 0 BOUNDARIES REGION 1 start ( -10, 0) natural(V) = 0 line to (10,0) value(V) = 0 line to (10,10) to (-10,10) to close REGION 2 start 'cw' (-3*rad,0) value(V)=-10 arc(center=-2*rad,0) to (-2*rad,rad) arc(center=-2*rad,0) to (-rad,0) natural(V)=0 line to close REGION 3 start 'ccw' (3*rad,0) value(V)=-10 arc(center=2*rad,0) to (2*rad,rad) arc(center=2*rad,0) to (rad,0) natural(V)=0 line to close REGION 4 start 'cw_el' (-6*rad,rad/100) value(V)=-10 arc(center=-5*rad,0) to (-5*rad,rad) arc(center=-5*rad,0) to (-4*rad,rad/100) natural(V)=0 line to close PLOTS contour(V) zoom(-6,0, 9,3) SUMMARY report( SURF_INTEGRAL(1,"cw")) as 'Clockwise Area (mm2)' report( SURF_INTEGRAL(1,"ccw")) as 'Counter-clockwise Area (mm2)' report( SURF_INTEGRAL(1,"cw_el")) as 'Clockwise, elevated Area (mm2)' END | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 738 Registered: 06-2003 |
FlexPDE appears to be counting bounding cells wrong, but I don't yet know the reason. I will continue to look into it. | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 741 Registered: 06-2003 |
The process by which this happens is a little convoluted, but it is clearly a hole in the logic FlexPDE applies to finding boundary integrals. Consider the case of a single ball with the boundary traced either counter-clockwise or clockwise (the case that produces the doubled integral). Then we have region 1 outside the ball and region 2 inside the ball. Your path walks the entire perimeter of the half circle, including the diameter. The diameter makes no contribution to the integral, since r=0 there. But.... In order to support the integration of functions that are discontinuous at boundaries, the boundary integral process needs to know what regions to evaluate the integrands in. By default, FlexPDE integrates in the material on the left side of the boundary as it is traversed. If the left side is exterior or if the right side region is specified for the evaluation, then the right side is used. In your case, the counterclockwise traversal has region 2 on its left for all segments of the boundary, so region 2 is the only region tagged for evaluation. Clockwise traversal selects region 1 for the arcs, but the diameter has no left side material, so it selects region 2 for evaluation. Unfortunately, these region selections are not tabulated by boundary segment, but only as a group for the integral path. So both region 1 and region 2 are selected for evaluation along the arcs as well as the diameter. Since there are two cells adjoining the integration path on the arcs, the integral is evaluated twice. This is a hole in the logic used by FlexPDE. I will fix it as soon as possible. But in the meantime, you can work around the problem by putting your names on Features that extend only over the part of the path that you want to integrate. In other words, remove the name from the boundary path and add a feature: Feature start 'cw' (-3*rad,0) arc(center=-2*rad,0) to (-2*rad,rad) to (-rad,0) This leaves the diameter out of the path and selects only region 1 for evaluation. The other canonical method for curing this, namely the explicit selection of evaluation region in the integral: surf_integral(1,"cw",2) also has logic troubles and can't find any cells to integrate in. I will fix both of these errors. | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 742 Registered: 06-2003 |
I have corrected the errors in mapping integration cells and posted a revised version for windows/32 at www.pdesolutions.com/download/xfpde5014x3win.exe This is a self-extracting archive containing only flexpde5.exe. It is not a full install. Extract over your old file, or change the name of the old file first. | ||
| morton777 New member Username: morton777 Post Number: 1 Registered: 10-2010 |
When calculating surface integrals, can I use the vector identity: |r_u X r_v| = |r_u|*|r_v|sin(theta) ? Also, if this is valid, what is the value of sin(theta)? | ||
| rgnelson Moderator Username: rgnelson Post Number: 1418 Registered: 06-2003 |
sin_theta = cross(r_u, r_v)/(magnitude(r_u)*magnitude(r_v)) |