[Paraview] Precisions about Surface Flow Filter

[Samuel Key]

All,

Users frequently want or call into question how surface facet areas are 
calculated.

For what it is worth, I have attached a closed-form (exact) calculation 
for bilinear 4-node quadrilateral facets in three-dimensions. The code 
'scrap' that implements this follows:
------------------------------------------------------------------
!!
!! Exact middle-surface area calculation: dA = (A,(P,Xi1 x P,Xi2))
!!
       Area = ZERO
       AK = SECTION_2D(SecID)%AK

       DO I = 1,4
         DO J = 1,4
           DO K = 1,4

             Earr = Ar(I)*Ps(J)*Pt(K) + As(I)*Pt(J)*Pr(K) + 
At(I)*Pr(J)*Ps(K)
&           - Ar(I)*Pt(J)*Ps(K) - As(I)*Pr(J)*Pt(K) - At(I)*Ps(J)*Pr(K)

             Area = Area + Earr*AK(I,J,K)
           ENDDO
         ENDDO
       ENDDO
------------------------------------------------------------------

Sam Key

On 2/11/2011 10:34 AM, Berk Geveci wrote:
> The easiest way to answer this is to show the code. Below is what it
> does for triangles. You can find what it does for other polygon types
> in here:
>
> http://paraview.org/gitweb?p=ParaView.git;a=blob_plain;f=Servers/Filters/vtkIntegrateAttributes.cxx;hb=HEAD
>
> Note that the actual filter is vtkIntegrateFlowThroughSurface which
> uses vtkSurfaceVectors to calculate vector . normal and then
> vtkIntegrateAttributes to integrate it over a polygonal surface. If
> you find a bug or inaccuracy that can be improved, please let us know.
>
> void vtkIntegrateAttributes::IntegrateTriangle(vtkDataSet* input,
>                                                 vtkUnstructuredGrid* output,
>                                                 vtkIdType cellId,
>                                                 vtkIdType pt1Id,
>                                                 vtkIdType pt2Id,
>                                                 vtkIdType pt3Id)
> {
>    double pt1[3], pt2[3], pt3[3];
>    double mid[3], v1[3], v2[3];
>    double cross[3];
>    double k;
>
>    input->GetPoint(pt1Id,pt1);
>    input->GetPoint(pt2Id,pt2);
>    input->GetPoint(pt3Id,pt3);
>
>    // Compute two legs.
>    v1[0] = pt2[0] - pt1[0];
>    v1[1] = pt2[1] - pt1[1];
>    v1[2] = pt2[2] - pt1[2];
>    v2[0] = pt3[0] - pt1[0];
>    v2[1] = pt3[1] - pt1[1];
>    v2[2] = pt3[2] - pt1[2];
>
>    // Use the cross product to compute the area of the parallelogram.
>    vtkMath::Cross(v1,v2,cross);
>    k = sqrt(cross[0]*cross[0] + cross[1]*cross[1] + cross[2]*cross[2]) * 0.5;
>
>    if (k == 0.0)
>      {
>      return;
>      }
>    this->Sum += k;
>
>    // Compute the middle, which is really just another attribute.
>    mid[0] = (pt1[0]+pt2[0]+pt3[0])/3.0;
>    mid[1] = (pt1[1]+pt2[1]+pt3[1])/3.0;
>    mid[2] = (pt1[2]+pt2[2]+pt3[2])/3.0;
>    // Add weighted to sumCenter.
>    this->SumCenter[0] += mid[0]*k;
>    this->SumCenter[1] += mid[1]*k;
>    this->SumCenter[2] += mid[2]*k;
>
>    // Now integrate the rest of the attributes.
>    this->IntegrateData3(input->GetPointData(), output->GetPointData(),
>                         pt1Id, pt2Id, pt3Id, k,
>                         *this->PointFieldList, this->FieldListIndex);
>    this->IntegrateData1(input->GetCellData(), output->GetCellData(), cellId, k,
>      *this->CellFieldList, this->FieldListIndex);
> }
>
> void vtkIntegrateAttributes::IntegrateData3(vtkDataSetAttributes* inda,
>    vtkDataSetAttributes* outda,
>    vtkIdType pt1Id, vtkIdType pt2Id,
>    vtkIdType pt3Id, double k,
>    vtkIntegrateAttributes::vtkFieldList&  fieldList, int index)
> {
>    int numArrays, i, numComponents, j;
>    vtkDataArray* inArray;
>    vtkDataArray* outArray;
>    numArrays = fieldList.GetNumberOfFields();
>    double vIn1, vIn2, vIn3, dv, vOut;
>    for (i = 0; i<  numArrays; ++i)
>      {
>      if (fieldList.GetFieldIndex(i)<  0)
>        {
>        continue;
>        }
>      // We could template for speed.
>      inArray = inda->GetArray(fieldList.GetDSAIndex(index, i));
>      outArray = outda->GetArray(fieldList.GetFieldIndex(i));
>      numComponents = inArray->GetNumberOfComponents();
>      for (j = 0; j<  numComponents; ++j)
>        {
>        vIn1 = inArray->GetComponent(pt1Id, j);
>        vIn2 = inArray->GetComponent(pt2Id, j);
>        vIn3 = inArray->GetComponent(pt3Id, j);
>        vOut = outArray->GetComponent(0, j);
>        dv = (vIn1+vIn2+vIn3)/3.0;
>        vOut += dv*k;
>        outArray->SetComponent(0, j, vOut);
>        }
>      }
> }
>
>
> 2011/2/11 Aurélien Marsan<aur.marsan at gmail.com>:
>> Dear all,
>>
>> Using the surface flow filter, I get results that are not those expected.
>> Could you give me some informations about how the integration is performed ?
>>
>> Is this a "local" approach, computing the surface of each face, and doing
>> some point data to face data interpolation ?
>>
>> Thanks,
>>
>> Aurélien.
>> _______________________________________________
>> Powered by www.kitware.com
>>
>> Visit other Kitware open-source projects at
>> http://www.kitware.com/opensource/opensource.html
>>
>> Please keep messages on-topic and check the ParaView Wiki at:
>> http://paraview.org/Wiki/ParaView
>>
>> Follow this link to subscribe/unsubscribe:
>> http://www.paraview.org/mailman/listinfo/paraview
>>
>>
> _______________________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html
>
> Please keep messages on-topic and check the ParaView Wiki at: http://paraview.org/Wiki/ParaView
>
> Follow this link to subscribe/unsubscribe:
> http://www.paraview.org/mailman/listinfo/paraview
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: el-area.pdf
Type: application/pdf
Size: 76744 bytes
Desc: not available
URL: <http://www.paraview.org/pipermail/paraview/attachments/20110211/58ba1141/attachment-0001.pdf>