From c150e1e4c40479c83b47cfb732f240c5df3b5d62 Mon Sep 17 00:00:00 2001 From: eap Date: Tue, 11 Apr 2017 19:40:35 +0300 Subject: [PATCH] 23422: EDF 14312 - Strange behavior of Viscous Layer Surface offset + smooth + Extract Morph from StdMeshers_Projection_2D.cxx + Fix formatting in SMESH_ScalarBarActor.cxx --- .../gui/SMESH/input/smeshpy_interface.doc | 4 +- src/OBJECT/SMESH_ScalarBarActor.cxx | 508 +++++++------- src/OBJECT/SMESH_ScalarBarActor.h | 2 +- src/SMESHGUI/SMESHGUI_DisplayEntitiesDlg.cxx | 2 +- src/StdMeshers/StdMeshers_ProjectionUtils.cxx | 320 ++++++++- src/StdMeshers/StdMeshers_ProjectionUtils.hxx | 32 +- src/StdMeshers/StdMeshers_Projection_2D.cxx | 238 +------ src/StdMeshers/StdMeshers_ViscousLayers.cxx | 660 +++++++++++++----- 8 files changed, 1101 insertions(+), 665 deletions(-) diff --git a/doc/salome/gui/SMESH/input/smeshpy_interface.doc b/doc/salome/gui/SMESH/input/smeshpy_interface.doc index 1dc47c407..524a4ceb0 100644 --- a/doc/salome/gui/SMESH/input/smeshpy_interface.doc +++ b/doc/salome/gui/SMESH/input/smeshpy_interface.doc @@ -2,10 +2,10 @@ \page smeshpy_interface_page Python interface -Python API for SALOME %Mesh module defines several classes that can +Python API of SALOME %Mesh module defines several classes that can be used for easy mesh creation and edition. -Documentation for SALOME %Mesh module Python API is available in two forms: +Documentation of SALOME %Mesh module Python API is available in two forms: - Structured documentation, where all methods and classes are grouped by their functionality. - Linear documentation grouped only by classes, declared diff --git a/src/OBJECT/SMESH_ScalarBarActor.cxx b/src/OBJECT/SMESH_ScalarBarActor.cxx index 1f5b00b63..c24369e0f 100644 --- a/src/OBJECT/SMESH_ScalarBarActor.cxx +++ b/src/OBJECT/SMESH_ScalarBarActor.cxx @@ -28,16 +28,16 @@ #include #include +#include #include #include #include +#include #include #include #include #include #include -#include -#include #define SHRINK_COEF 0.08; @@ -51,13 +51,14 @@ vtkCxxSetObjectMacro(SMESH_ScalarBarActor,TitleTextProperty,vtkTextProperty); // Instantiate object with 64 maximum colors; 5 labels; %%-#6.3g label // format, no title, and vertical orientation. The initial scalar bar // size is (0.05 x 0.8) of the viewport size. -SMESH_ScalarBarActor::SMESH_ScalarBarActor() { +SMESH_ScalarBarActor::SMESH_ScalarBarActor() +{ this->LookupTable = NULL; this->Position2Coordinate->SetValue(0.17, 0.8); - + this->PositionCoordinate->SetCoordinateSystemToNormalizedViewport(); this->PositionCoordinate->SetValue(0.82,0.1); - + this->MaximumNumberOfColors = 64; this->NumberOfLabels = 5; this->NumberOfLabelsBuilt = 0; @@ -74,7 +75,7 @@ SMESH_ScalarBarActor::SMESH_ScalarBarActor() { this->TitleTextProperty = vtkTextProperty::New(); this->TitleTextProperty->ShallowCopy(this->LabelTextProperty); - this->LabelFormat = new char[8]; + this->LabelFormat = new char[8]; sprintf(this->LabelFormat,"%s","%-#6.3g"); this->TitleMapper = vtkTextMapper::New(); @@ -82,7 +83,7 @@ SMESH_ScalarBarActor::SMESH_ScalarBarActor() { this->TitleActor->SetMapper(this->TitleMapper); this->TitleActor->GetPositionCoordinate()-> SetReferenceCoordinate(this->PositionCoordinate); - + this->TextMappers = NULL; this->TextActors = NULL; @@ -104,7 +105,7 @@ SMESH_ScalarBarActor::SMESH_ScalarBarActor() { myDistribution = vtkPolyData::New(); myDistributionMapper = vtkPolyDataMapper2D::New(); myDistributionMapper->SetInputData(this->myDistribution); - + myDistributionActor = vtkActor2D::New(); myDistributionActor->SetMapper(this->myDistributionMapper); myDistributionActor->GetPositionCoordinate()-> @@ -129,12 +130,12 @@ void SMESH_ScalarBarActor::ReleaseGraphicsResources(vtkWindow *win) { this->TitleActor->ReleaseGraphicsResources(win); if (this->TextMappers != NULL ) - { + { for (int i=0; i < this->NumberOfLabelsBuilt; i++) - { + { this->TextActors[i]->ReleaseGraphicsResources(win); - } } + } this->ScalarBarActor->ReleaseGraphicsResources(win); // rnv begin // Customization of the vtkScalarBarActor to show distribution histogram. @@ -143,41 +144,42 @@ void SMESH_ScalarBarActor::ReleaseGraphicsResources(vtkWindow *win) /*--------------------------------------------------------------------------*/ -SMESH_ScalarBarActor::~SMESH_ScalarBarActor() { - if (this->LabelFormat) - { +SMESH_ScalarBarActor::~SMESH_ScalarBarActor() +{ + if (this->LabelFormat) + { delete [] this->LabelFormat; this->LabelFormat = NULL; - } + } this->TitleMapper->Delete(); this->TitleActor->Delete(); if (this->TextMappers != NULL ) - { + { for (int i=0; i < this->NumberOfLabelsBuilt; i++) - { + { this->TextMappers[i]->Delete(); this->TextActors[i]->Delete(); - } + } delete [] this->TextMappers; delete [] this->TextActors; - } + } this->ScalarBar->Delete(); this->ScalarBarMapper->Delete(); this->ScalarBarActor->Delete(); if (this->Title) - { + { delete [] this->Title; this->Title = NULL; - } - + } + this->SetLookupTable(NULL); this->SetLabelTextProperty(NULL); this->SetTitleTextProperty(NULL); - + // rnv begin // Customization of the vtkScalarBarActor to show distribution histogram: myDistribution->Delete(); @@ -191,26 +193,26 @@ int SMESH_ScalarBarActor::RenderOverlay(vtkViewport *viewport) { int renderedSomething = 0; int i; - + // Everything is built, just have to render if (this->Title != NULL) - { + { renderedSomething += this->TitleActor->RenderOverlay(viewport); - } - if (!myTitleOnlyVisibility) { + } + if ( !myTitleOnlyVisibility ) { this->ScalarBarActor->RenderOverlay(viewport); this->myDistributionActor->RenderOverlay(viewport); - if( this->TextActors == NULL) - { - vtkWarningMacro(<<"Need a mapper to render a scalar bar"); - return renderedSomething; - } - - for (i=0; iNumberOfLabels; i++) - { + if ( this->TextActors == NULL ) + { + vtkWarningMacro(<<"Need a mapper to render a scalar bar"); + return renderedSomething; + } + + for ( i=0; iNumberOfLabels; i++ ) + { renderedSomething += this->TextActors[i]->RenderOverlay(viewport); - } - } + } + } renderedSomething = (renderedSomething > 0)?(1):(0); return renderedSomething; @@ -223,75 +225,75 @@ int SMESH_ScalarBarActor::RenderOpaqueGeometry(vtkViewport *viewport) int renderedSomething = 0; int i; int size[2]; - + if (!this->LookupTable) - { + { vtkWarningMacro(<<"Need a mapper to render a scalar bar"); return 0; - } + } if (!this->TitleTextProperty) - { + { vtkErrorMacro(<<"Need title text property to render a scalar bar"); return 0; - } + } if (!this->LabelTextProperty) - { + { vtkErrorMacro(<<"Need label text property to render a scalar bar"); return 0; - } + } // Check to see whether we have to rebuild everything int positionsHaveChanged = 0; - if (viewport->GetMTime() > this->BuildTime || - (viewport->GetVTKWindow() && + if (viewport->GetMTime() > this->BuildTime || + (viewport->GetVTKWindow() && viewport->GetVTKWindow()->GetMTime() > this->BuildTime)) - { + { // if the viewport has changed we may - or may not need // to rebuild, it depends on if the projected coords chage int *barOrigin; barOrigin = this->PositionCoordinate->GetComputedViewportValue(viewport); - size[0] = + size[0] = this->Position2Coordinate->GetComputedViewportValue(viewport)[0] - barOrigin[0]; - size[1] = + size[1] = this->Position2Coordinate->GetComputedViewportValue(viewport)[1] - barOrigin[1]; - if (this->LastSize[0] != size[0] || + if (this->LastSize[0] != size[0] || this->LastSize[1] != size[1] || - this->LastOrigin[0] != barOrigin[0] || + this->LastOrigin[0] != barOrigin[0] || this->LastOrigin[1] != barOrigin[1]) - { + { positionsHaveChanged = 1; - } } - + } + // Check to see whether we have to rebuild everything - if (positionsHaveChanged || - this->GetMTime() > this->BuildTime || - this->LookupTable->GetMTime() > this->BuildTime || - this->LabelTextProperty->GetMTime() > this->BuildTime || - this->TitleTextProperty->GetMTime() > this->BuildTime) - { + if ( positionsHaveChanged || + this->GetMTime() > this->BuildTime || + this->LookupTable->GetMTime() > this->BuildTime || + this->LabelTextProperty->GetMTime() > this->BuildTime || + this->TitleTextProperty->GetMTime() > this->BuildTime) + { vtkDebugMacro(<<"Rebuilding subobjects"); // Delete previously constructed objects // - if (this->TextMappers != NULL ) + if ( this->TextMappers != NULL ) + { + for ( i = 0; i < this->NumberOfLabelsBuilt; i++ ) { - for (i=0; i < this->NumberOfLabelsBuilt; i++) - { this->TextMappers[i]->Delete(); this->TextActors[i]->Delete(); - } + } delete [] this->TextMappers; delete [] this->TextActors; - } + } // Build scalar bar object; determine its type // - // is this a vtkLookupTable or a subclass of vtkLookupTable + // is this a vtkLookupTable or a subclass of vtkLookupTable // with its scale set to log // NOTE: it's possible we could to without the 'lut' variable // later in the code, but if the vtkLookupTableSafeDownCast operation @@ -300,13 +302,13 @@ int SMESH_ScalarBarActor::RenderOpaqueGeometry(vtkViewport *viewport) vtkLookupTable *LUT = vtkLookupTable::SafeDownCast( this->LookupTable ); int isLogTable = 0; if ( LUT ) - { + { if ( LUT->GetScale() == VTK_SCALE_LOG10 ) - { - isLogTable = 1; - } + { + isLogTable = 1; } - + } + // we hard code how many steps to display vtkScalarsToColors *lut = this->LookupTable; int numColors = this->MaximumNumberOfColors; @@ -332,7 +334,8 @@ int SMESH_ScalarBarActor::RenderOpaqueGeometry(vtkViewport *viewport) if(!distrVisibility) vtkDebugMacro(<<" Distribution invisible, because numColors == this->myNbValues.size()"); - if ( distrVisibility && GetDistributionVisibility() ) { + if ( distrVisibility && GetDistributionVisibility() ) + { for ( i = 0 ; i < (int)myNbValues.size(); i++ ) { if ( myNbValues[i] ) { numPositiveVal++; @@ -349,16 +352,21 @@ int SMESH_ScalarBarActor::RenderOpaqueGeometry(vtkViewport *viewport) this->myDistribution->SetPolys(distrPolys); distrPts->Delete(); distrPolys->Delete(); - if ( myDistributionColoringType == SMESH_MULTICOLOR_TYPE ) { + if ( myDistributionColoringType == SMESH_MULTICOLOR_TYPE ) + { distColors = vtkUnsignedCharArray::New(); distColors->SetNumberOfComponents(3); distColors->SetNumberOfTuples(numPositiveVal); this->myDistribution->GetCellData()->SetScalars(distColors); distColors->Delete(); - } else if( myDistributionColoringType == SMESH_MONOCOLOR_TYPE ){ + } + else if( myDistributionColoringType == SMESH_MONOCOLOR_TYPE ) + { this->myDistribution->GetCellData()->SetScalars(NULL); } - } else { + } + else + { myDistribution->Reset(); } // rnv end @@ -373,22 +381,22 @@ int SMESH_ScalarBarActor::RenderOpaqueGeometry(vtkViewport *viewport) // get the viewport size in display coordinates int *barOrigin, barWidth, barHeight, distrHeight; barOrigin = this->PositionCoordinate->GetComputedViewportValue(viewport); - size[0] = + size[0] = this->Position2Coordinate->GetComputedViewportValue(viewport)[0] - barOrigin[0]; - size[1] = + size[1] = this->Position2Coordinate->GetComputedViewportValue(viewport)[1] - barOrigin[1]; this->LastOrigin[0] = barOrigin[0]; this->LastOrigin[1] = barOrigin[1]; this->LastSize[0] = size[0]; this->LastSize[1] = size[1]; - + // Update all the composing objects this->TitleActor->SetProperty(this->GetProperty()); this->TitleMapper->SetInput(this->Title); if (this->TitleTextProperty->GetMTime() > this->BuildTime) - { + { // Shallow copy here so that the size of the title prop is not affected // by the automatic adjustment of its text mapper's size (i.e. its // mapper's text property is identical except for the font size @@ -397,17 +405,17 @@ int SMESH_ScalarBarActor::RenderOpaqueGeometry(vtkViewport *viewport) // the title and label text prop to be the same. this->TitleMapper->GetTextProperty()->ShallowCopy(this->TitleTextProperty); this->TitleMapper->GetTextProperty()->SetJustificationToCentered(); - } - + } + // find the best size for the title font int titleSize[2]; this->SizeTitle(titleSize, size, viewport); - + // find the best size for the ticks int labelSize[2]; this->AllocateAndSizeLabels(labelSize, size, viewport,range); this->NumberOfLabelsBuilt = this->NumberOfLabels; - + // generate points double x[3]; x[2] = 0.0; double delta, itemH, shrink; @@ -426,7 +434,7 @@ int SMESH_ScalarBarActor::RenderOpaqueGeometry(vtkViewport *viewport) barHeight = (int)(0.86*size[1]); delta=(double)barHeight/numColors; - + for ( i=0; iSetPoint(2*i+1,x); } - if(GetDistributionVisibility() && distrVisibility) { - // Distribution points + if ( GetDistributionVisibility() && distrVisibility ) { + // Distribution points shrink = delta*SHRINK_COEF; vtkIdType distPtsId=0; vtkIdType distPtsIds[4]; - for(i=0; iSetPoint(distPtsId++,x); @@ -462,7 +470,7 @@ int SMESH_ScalarBarActor::RenderOpaqueGeometry(vtkViewport *viewport) x[1] = i*delta+delta-shrink; // third point of polygon (quadrangle) - x[0] = 0; + x[0] = 0; distPtsIds[3] = distPtsId; distrPts->SetPoint(distPtsId++,x); @@ -475,19 +483,19 @@ int SMESH_ScalarBarActor::RenderOpaqueGeometry(vtkViewport *viewport) distrPolys->InsertNextCell(4,distPtsIds); } } - } + } } // rnv end else { barWidth = size[0]; - + // rnv begin // Customization of the vtkScalarBarActor to show distribution histogram. double coef1, delimeter=0.0; - if(GetDistributionVisibility() && distrVisibility) { + if ( GetDistributionVisibility() && distrVisibility ) { coef1=0.62; distrHeight = (int)((coef1/2)*size[1]); - //delimeter between distribution diagram and scalar bar + //delimeter between distribution diagram and scalar bar delimeter=0.02*size[1]; } else { @@ -495,64 +503,64 @@ int SMESH_ScalarBarActor::RenderOpaqueGeometry(vtkViewport *viewport) barHeight = (int)(coef1*size[1]); distrHeight = 0; } - + barHeight = (int)(coef1*size[1]); - + delta=(double)barWidth/numColors; - for (i=0; iSetPoint(2*i,x); + pts->SetPoint(2*i,x); x[1] = distrHeight + delimeter; pts->SetPoint(2*i+1,x); } - - if(GetDistributionVisibility() && distrVisibility) { - // Distribution points + + if ( GetDistributionVisibility() && distrVisibility ) { + // Distribution points shrink = delta*SHRINK_COEF; vtkIdType distPtsId=0; vtkIdType distPtsIds[4]; - for(i=0; iSetPoint(distPtsId++,x); - + // second point of polygon (quadrangle) - x[0] = i*delta+shrink; + x[0] = i*delta+shrink; x[1] = itemH; distPtsIds[3] = distPtsId; distrPts->SetPoint(distPtsId++,x); - + // third point of polygon (quadrangle) - x[0] = i*delta+delta-shrink; + x[0] = i*delta+delta-shrink; x[1] = 0; distPtsIds[1] = distPtsId; distrPts->SetPoint(distPtsId++,x); - + // fourth point of polygon (quadrangle) - x[0] = i*delta+delta-shrink; + x[0] = i*delta+delta-shrink; x[1] = itemH; distPtsIds[2] = distPtsId; distrPts->SetPoint(distPtsId++,x); - + // Add polygon into poly data distrPolys->InsertNextCell(4,distPtsIds); } - } + } } // rnv end } - + //polygons & cell colors unsigned char *rgba, *rgb; vtkIdType ptIds[4], dcCount=0; - for (i=0; iInsertNextCell(4,ptIds); if ( isLogTable ) - { - double rgbval = log10(range[0]) + + { + double rgbval = log10(range[0]) + i*(log10(range[1])-log10(range[0]))/(numColors -1); rgba = lut->MapValue(pow(10.0,rgbval)); - } + } else - { + { rgba = lut->MapValue(range[0] + (range[1] - range[0])* ((double)i /(numColors-1.0))); - } + } rgb = colors->GetPointer(3*i); //write into array directly rgb[0] = rgba[0]; rgb[1] = rgba[1]; rgb[2] = rgba[2]; - + // rnv begin // Customization of the vtkScalarBarActor to show distribution histogram. - if(myDistributionColoringType == SMESH_MULTICOLOR_TYPE && GetDistributionVisibility() && distrVisibility) - { - rgb = distColors->GetPointer(3*dcCount); //write into array directly - rgb[0] = rgba[0]; - rgb[1] = rgba[1]; - rgb[2] = rgba[2]; - dcCount++; - } + if ( myDistributionColoringType == SMESH_MULTICOLOR_TYPE && + GetDistributionVisibility() && + distrVisibility ) + { + rgb = distColors->GetPointer(3*dcCount); //write into array directly + rgb[0] = rgba[0]; + rgb[1] = rgba[1]; + rgb[2] = rgba[2]; + dcCount++; } + } // Now position everything properly // double val; - if (this->Orientation == VTK_ORIENT_VERTICAL) - { + if ( this->Orientation == VTK_ORIENT_VERTICAL ) + { int sizeTextData[2]; - + // center the title this->TitleActor->SetPosition(size[0]/2, 0.9*size[1]); - - for (i=0; i < this->NumberOfLabels; i++) + + for ( i = 0; i < this->NumberOfLabels; i++ ) + { + if ( this->NumberOfLabels > 1 ) { - if (this->NumberOfLabels > 1) - { val = (double)i/(this->NumberOfLabels-1) *barHeight; - } - else - { + } + else + { val = 0.5*barHeight; - } + } this->TextMappers[i]->GetSize(viewport,sizeTextData); this->TextMappers[i]->GetTextProperty()->SetJustificationToLeft(); this->TextActors[i]->SetPosition(barWidth+3, val - sizeTextData[1]/2); - } } + } else - { - this->TitleActor->SetPosition(size[0]/2, + { + this->TitleActor->SetPosition(size[0]/2, barHeight + labelSize[1] + 0.1*size[1]); - for (i=0; i < this->NumberOfLabels; i++) - { + for ( i = 0; i < this->NumberOfLabels; i++ ) + { this->TextMappers[i]->GetTextProperty()->SetJustificationToCentered(); if (this->NumberOfLabels > 1) - { + { val = (double)i/(this->NumberOfLabels-1) * barWidth; - } + } else - { + { val = 0.5*barWidth; - } - this->TextActors[i]->SetPosition(val, barHeight + 0.05*size[1]); } + this->TextActors[i]->SetPosition(val, barHeight + 0.05*size[1]); } + } this->BuildTime.Modified(); - } + } // Everything is built, just have to render - if (this->Title != NULL) - { + if ( this->Title != NULL ) + { renderedSomething += this->TitleActor->RenderOpaqueGeometry(viewport); - } + } this->ScalarBarActor->RenderOpaqueGeometry(viewport); this->myDistributionActor->RenderOpaqueGeometry(viewport); - for (i=0; iNumberOfLabels; i++) - { + for ( i = 0; i < this->NumberOfLabels; i++ ) + { renderedSomething += this->TextActors[i]->RenderOpaqueGeometry(viewport); - } + } renderedSomething = (renderedSomething > 0)?(1):(0); @@ -659,50 +669,50 @@ void SMESH_ScalarBarActor::PrintSelf(ostream& os, vtkIndent indent) this->Superclass::PrintSelf(os,indent); if ( this->LookupTable ) - { + { os << indent << "Lookup Table:\n"; this->LookupTable->PrintSelf(os,indent.GetNextIndent()); - } + } else - { + { os << indent << "Lookup Table: (none)\n"; - } + } if (this->TitleTextProperty) - { + { os << indent << "Title Text Property:\n"; this->TitleTextProperty->PrintSelf(os,indent.GetNextIndent()); - } + } else - { + { os << indent << "Title Text Property: (none)\n"; - } + } if (this->LabelTextProperty) - { + { os << indent << "Label Text Property:\n"; this->LabelTextProperty->PrintSelf(os,indent.GetNextIndent()); - } + } else - { + { os << indent << "Label Text Property: (none)\n"; - } + } os << indent << "Title: " << (this->Title ? this->Title : "(none)") << "\n"; - os << indent << "Maximum Number Of Colors: " + os << indent << "Maximum Number Of Colors: " << this->MaximumNumberOfColors << "\n"; os << indent << "Number Of Labels: " << this->NumberOfLabels << "\n"; os << indent << "Number Of Labels Built: " << this->NumberOfLabelsBuilt << "\n"; os << indent << "Orientation: "; if ( this->Orientation == VTK_ORIENT_HORIZONTAL ) - { + { os << "Horizontal\n"; - } + } else - { + { os << "Vertical\n"; - } + } os << indent << "Label Format: " << this->LabelFormat << "\n"; } @@ -712,7 +722,7 @@ void SMESH_ScalarBarActor::ShallowCopy(vtkProp *prop) { SMESH_ScalarBarActor *a = SMESH_ScalarBarActor::SafeDownCast(prop); if ( a != NULL ) - { + { this->SetPosition2(a->GetPosition2()); this->SetLookupTable(a->GetLookupTable()); this->SetMaximumNumberOfColors(a->GetMaximumNumberOfColors()); @@ -721,25 +731,25 @@ void SMESH_ScalarBarActor::ShallowCopy(vtkProp *prop) this->SetTitleTextProperty(a->GetTitleTextProperty()); this->SetLabelFormat(a->GetLabelFormat()); this->SetTitle(a->GetTitle()); - this->GetPositionCoordinate()->SetCoordinateSystem( - a->GetPositionCoordinate()->GetCoordinateSystem()); - this->GetPositionCoordinate()->SetValue( - a->GetPositionCoordinate()->GetValue()); - this->GetPosition2Coordinate()->SetCoordinateSystem( - a->GetPosition2Coordinate()->GetCoordinateSystem()); - this->GetPosition2Coordinate()->SetValue( - a->GetPosition2Coordinate()->GetValue()); - } + this->GetPositionCoordinate()->SetCoordinateSystem + (a->GetPositionCoordinate()->GetCoordinateSystem()); + this->GetPositionCoordinate()->SetValue + (a->GetPositionCoordinate()->GetValue()); + this->GetPosition2Coordinate()->SetCoordinateSystem + (a->GetPosition2Coordinate()->GetCoordinateSystem()); + this->GetPosition2Coordinate()->SetValue + (a->GetPosition2Coordinate()->GetValue()); + } // Now do superclass this->vtkActor2D::ShallowCopy(prop); } //---------------------------------------------------------------------------- -void SMESH_ScalarBarActor::AllocateAndSizeLabels(int *labelSize, - int *size, - vtkViewport *viewport, - double *range) +void SMESH_ScalarBarActor::AllocateAndSizeLabels(int *labelSize, + int *size, + vtkViewport *viewport, + double *range) { labelSize[0] = labelSize[1] = 0; @@ -750,55 +760,55 @@ void SMESH_ScalarBarActor::AllocateAndSizeLabels(int *labelSize, double val; int i; - + // TODO: this should be optimized, maybe by keeping a list of // allocated mappers, in order to avoid creation/destruction of // their underlying text properties (i.e. each time a mapper is // created, text properties are created and shallow-assigned a font size // which value might be "far" from the target font size). - // is this a vtkLookupTable or a subclass of vtkLookupTable + // is this a vtkLookupTable or a subclass of vtkLookupTable // with its scale set to log vtkLookupTable *LUT = vtkLookupTable::SafeDownCast( this->LookupTable ); int isLogTable = 0; if ( LUT ) - { + { if ( LUT->GetScale() == VTK_SCALE_LOG10 ) - { - isLogTable = 1; - } + { + isLogTable = 1; } + } - for (i=0; i < this->NumberOfLabels; i++) - { + for ( i = 0; i < this->NumberOfLabels; i++ ) + { this->TextMappers[i] = vtkTextMapper::New(); if ( isLogTable ) - { + { double lval; - if (this->NumberOfLabels > 1) - { + if ( this->NumberOfLabels > 1 ) + { lval = log10(range[0]) + (double)i/(this->NumberOfLabels-1) * (log10(range[1])-log10(range[0])); - } + } else - { + { lval = log10(range[0]) + 0.5*(log10(range[1])-log10(range[0])); - } - val = pow(10.0,lval); } + val = pow(10.0,lval); + } else + { + if ( this->NumberOfLabels > 1 ) { - if (this->NumberOfLabels > 1) - { - val = range[0] + + val = range[0] + (double)i/(this->NumberOfLabels-1) * (range[1]-range[0]); - } + } else - { + { val = range[0] + 0.5*(range[1]-range[0]); - } } + } sprintf(string, this->LabelFormat, val); this->TextMappers[i]->SetInput(string); @@ -809,64 +819,63 @@ void SMESH_ScalarBarActor::AllocateAndSizeLabels(int *labelSize, // which will be modified later). This allows text actors to // share the same text property, and in that case specifically allows // the title and label text prop to be the same. - this->TextMappers[i]->GetTextProperty()->ShallowCopy( - this->LabelTextProperty); + this->TextMappers[i]->GetTextProperty()->ShallowCopy(this->LabelTextProperty); this->TextActors[i] = vtkActor2D::New(); this->TextActors[i]->SetMapper(this->TextMappers[i]); this->TextActors[i]->SetProperty(this->GetProperty()); this->TextActors[i]->GetPositionCoordinate()-> SetReferenceCoordinate(this->PositionCoordinate); - } + } - if (this->NumberOfLabels) - { + if ( this->NumberOfLabels ) + { int targetWidth, targetHeight; // rnv begin // Customization of the vtkScalarBarActor to show distribution histogram. bool distrVisibility = ( this->MaximumNumberOfColors == (int) this->myNbValues.size() ); double coef; - if( GetDistributionVisibility() && distrVisibility ) - if(this->Orientation == VTK_ORIENT_VERTICAL) + if ( GetDistributionVisibility() && distrVisibility ) + if ( this->Orientation == VTK_ORIENT_VERTICAL ) coef = 0.4; - else + else coef = 0.18; - else - if(this->Orientation == VTK_ORIENT_VERTICAL) + else + if (this->Orientation == VTK_ORIENT_VERTICAL ) coef = 0.6; - else + else coef=0.25; if ( this->Orientation == VTK_ORIENT_VERTICAL ) - { + { targetWidth = (int)(coef*size[0]); targetHeight = (int)(0.86*size[1]/this->NumberOfLabels); - } + } else - { + { targetWidth = (int)(size[0]*0.8/this->NumberOfLabels); targetHeight = (int)(coef*size[1]); - } - // rnv end - - vtkTextMapper::SetMultipleConstrainedFontSize(viewport, - targetWidth, - targetHeight, - this->TextMappers, - this->NumberOfLabels, - labelSize); } + // rnv end + + vtkTextMapper::SetMultipleConstrainedFontSize( viewport, + targetWidth, + targetHeight, + this->TextMappers, + this->NumberOfLabels, + labelSize ); + } } //---------------------------------------------------------------------------- -void SMESH_ScalarBarActor::SizeTitle(int *titleSize, - int *size, +void SMESH_ScalarBarActor::SizeTitle(int *titleSize, + int *size, vtkViewport *viewport) { titleSize[0] = titleSize[1] = 0; - if (this->Title == NULL || !strlen(this->Title)) + if ( this->Title == NULL || !strlen(this->Title) ) { return; } @@ -899,36 +908,43 @@ void SMESH_ScalarBarActor::SizeTitle(int *titleSize, /*--------------------------------------------------------------------------*/ -void SMESH_ScalarBarActor::SetDistributionVisibility(int flag) { - myDistributionActor->SetVisibility(flag); +void SMESH_ScalarBarActor::SetDistributionVisibility( int flag ) +{ + myDistributionActor->SetVisibility( flag ); Modified(); } /*--------------------------------------------------------------------------*/ -int SMESH_ScalarBarActor::GetDistributionVisibility() { +int SMESH_ScalarBarActor::GetDistributionVisibility() +{ return myDistributionActor->GetVisibility(); } -void SMESH_ScalarBarActor::SetDistribution(std::vector theNbValues) { +void SMESH_ScalarBarActor::SetDistribution( const std::vector& theNbValues ) +{ myNbValues = theNbValues; -} +} -void SMESH_ScalarBarActor::SetDistributionColor (double rgb[3]) { +void SMESH_ScalarBarActor::SetDistributionColor( double rgb[3] ) +{ myDistributionActor->GetProperty()->SetColor(rgb); Modified(); } -void SMESH_ScalarBarActor::GetDistributionColor (double rgb[3]) { +void SMESH_ScalarBarActor::GetDistributionColor( double rgb[3] ) +{ myDistributionActor->GetProperty()->GetColor(rgb); } -void SMESH_ScalarBarActor::SetTitleOnlyVisibility( bool theTitleOnlyVisibility) { +void SMESH_ScalarBarActor::SetTitleOnlyVisibility( bool theTitleOnlyVisibility ) +{ myTitleOnlyVisibility = theTitleOnlyVisibility; } -bool SMESH_ScalarBarActor::GetTitleOnlyVisibility() { +bool SMESH_ScalarBarActor::GetTitleOnlyVisibility() +{ return myTitleOnlyVisibility; } diff --git a/src/OBJECT/SMESH_ScalarBarActor.h b/src/OBJECT/SMESH_ScalarBarActor.h index ffc91bb11..c5e721b23 100644 --- a/src/OBJECT/SMESH_ScalarBarActor.h +++ b/src/OBJECT/SMESH_ScalarBarActor.h @@ -174,7 +174,7 @@ class SMESHOBJECT_EXPORT SMESH_ScalarBarActor: public vtkActor2D { virtual int GetDistributionVisibility(); // Description: // Set distribution - virtual void SetDistribution(std::vector theNbValues); + virtual void SetDistribution(const std::vector& theNbValues); // Description: // Set distribution coloring type (SMESH_MONOCOLOR_TYPE or SMESH_MULTICOLOR_TYPE) diff --git a/src/SMESHGUI/SMESHGUI_DisplayEntitiesDlg.cxx b/src/SMESHGUI/SMESHGUI_DisplayEntitiesDlg.cxx index 5af7f7815..14088fe1a 100644 --- a/src/SMESHGUI/SMESHGUI_DisplayEntitiesDlg.cxx +++ b/src/SMESHGUI/SMESHGUI_DisplayEntitiesDlg.cxx @@ -97,7 +97,7 @@ SMESHGUI_DisplayEntitiesDlg::SMESHGUI_DisplayEntitiesDlg( QWidget* parent ) hl->addWidget( nb0DElemsLab, 0, 1 ); my0DElemsTB->setEnabled( nbElements ); nb0DElemsLab->setEnabled( nbElements ); - myNbTypes += ( nbElements > 0 ); + myNbTypes = ( nbElements > 0 ); // Edges nbElements = myActor ? myActor->GetObject()->GetNbEntities( SMDSAbs_Edge ) : aMesh->NbEdges(); diff --git a/src/StdMeshers/StdMeshers_ProjectionUtils.cxx b/src/StdMeshers/StdMeshers_ProjectionUtils.cxx index a0164ef01..e7051db9b 100644 --- a/src/StdMeshers/StdMeshers_ProjectionUtils.cxx +++ b/src/StdMeshers/StdMeshers_ProjectionUtils.cxx @@ -28,6 +28,7 @@ #include "StdMeshers_ProjectionUtils.hxx" #include "SMDS_EdgePosition.hxx" +#include "SMDS_FacePosition.hxx" #include "SMESHDS_Mesh.hxx" #include "SMESH_Algo.hxx" #include "SMESH_Block.hxx" @@ -46,6 +47,7 @@ #include "utilities.h" #include +#include #include #include #include @@ -2787,4 +2789,320 @@ namespace StdMeshers_ProjectionUtils } return true; } -} + + //================================================================================ + /*! + * \brief Add in-FACE nodes surrounding a given node to a queue + */ + //================================================================================ + + typedef list< pair< const SMDS_MeshNode*, const BRepMesh_Triangle* > > TNodeTriaList; + + void addCloseNodes( const SMDS_MeshNode* srcNode, + const BRepMesh_Triangle* bmTria, + const int srcFaceID, + TNodeTriaList & noTriQueue ) + { + // find in-FACE nodes + SMDS_ElemIteratorPtr elems = srcNode->GetInverseElementIterator(SMDSAbs_Face); + while ( elems->more() ) + { + const SMDS_MeshElement* elem = elems->next(); + if ( elem->getshapeId() == srcFaceID ) + { + for ( int i = 0, nb = elem->NbNodes(); i < nb; ++i ) + { + const SMDS_MeshNode* n = elem->GetNode( i ); + if ( !n->isMarked() ) + noTriQueue.push_back( make_pair( n, bmTria )); + } + } + } + } + + //================================================================================ + /*! + * \brief Find a delauney triangle containing a given 2D point and return + * barycentric coordinates within the found triangle + */ + //================================================================================ + + const BRepMesh_Triangle* findTriangle( const gp_XY& uv, + const BRepMesh_Triangle* bmTria, + Handle(BRepMesh_DataStructureOfDelaun)& triaDS, + double bc[3] ) + { + int nodeIDs[3]; + gp_XY nodeUVs[3]; + int linkIDs[3]; + Standard_Boolean ori[3]; + + while ( bmTria ) + { + // check bmTria + + triaDS->ElementNodes( *bmTria, nodeIDs ); + nodeUVs[0] = triaDS->GetNode( nodeIDs[0] ).Coord(); + nodeUVs[1] = triaDS->GetNode( nodeIDs[1] ).Coord(); + nodeUVs[2] = triaDS->GetNode( nodeIDs[2] ).Coord(); + + SMESH_MeshAlgos::GetBarycentricCoords( uv, + nodeUVs[0], nodeUVs[1], nodeUVs[2], + bc[0], bc[1] ); + if ( bc[0] >= 0 && bc[1] >= 0 && bc[0] + bc[1] <= 1 ) + { + bc[2] = 1 - bc[0] - bc[1]; + return bmTria; + } + + // look for a neighbor triangle, which is adjacent to a link intersected + // by a segment( triangle center -> uv ) + + gp_XY gc = ( nodeUVs[0] + nodeUVs[1] + nodeUVs[2] ) / 3.; + gp_XY seg = uv - gc; + + bmTria->Edges( linkIDs, ori ); + int triaID = triaDS->IndexOf( *bmTria ); + bmTria = 0; + + for ( int i = 0; i < 3; ++i ) + { + const BRepMesh_PairOfIndex & triIDs = triaDS->ElementsConnectedTo( linkIDs[i] ); + if ( triIDs.Extent() < 2 ) + continue; // no neighbor triangle + + // check if a link intersects gc2uv + const BRepMesh_Edge & link = triaDS->GetLink( linkIDs[i] ); + const BRepMesh_Vertex & n1 = triaDS->GetNode( link.FirstNode() ); + const BRepMesh_Vertex & n2 = triaDS->GetNode( link.LastNode() ); + gp_XY uv1 = n1.Coord(); + gp_XY lin = n2.Coord() - uv1; // link direction + + double crossSegLin = seg ^ lin; + if ( Abs( crossSegLin ) < std::numeric_limits::min() ) + continue; // parallel + + double uSeg = ( uv1 - gc ) ^ lin / crossSegLin; + if ( 0. <= uSeg && uSeg <= 1. ) + { + bmTria = & triaDS->GetElement( triIDs.Index( 1 + ( triIDs.Index(1) == triaID ))); + break; + } + } + } + return bmTria; + } + + //================================================================================ + /*! + * \brief triangulate the srcFace in 2D + * \param [in] srcWires - boundary of the src FACE + */ + //================================================================================ + + Morph::Morph(const TSideVector& srcWires) + { + _srcSubMesh = srcWires[0]->GetMesh()->GetSubMesh( srcWires[0]->Face() ); + + // compute _scale + { + BRepAdaptor_Surface surf( srcWires[0]->Face() ); + const int nbDiv = 100; + const double uRange = surf.LastUParameter() - surf.FirstUParameter(); + const double vRange = surf.LastVParameter() - surf.FirstVParameter(); + const double dU = uRange / nbDiv; + const double dV = vRange / nbDiv; + double u = surf.FirstUParameter(), v = surf.FirstVParameter(); + gp_Pnt p0U = surf.Value( u, v ), p0V = p0U; + double lenU = 0, lenV = 0; + for ( ; u < surf.LastUParameter(); u += dU, v += dV ) + { + gp_Pnt p1U = surf.Value( u, surf.FirstVParameter() ); + lenU += p1U.Distance( p0U ); + p0U = p1U; + gp_Pnt p1V = surf.Value( surf.FirstUParameter(), v ); + lenV += p1V.Distance( p0V ); + p0V = p1V; + } + _scale.SetCoord( lenU / uRange, lenV / vRange ); + } + + // count boundary points + int iP = 1, nbP = 0; + for ( size_t iW = 0; iW < srcWires.size(); ++iW ) + nbP += srcWires[iW]->NbPoints() - 1; // 1st and last points coincide + + _bndSrcNodes.resize( nbP + 1 ); _bndSrcNodes[0] = 0; + + // fill boundary points + BRepMesh::Array1OfVertexOfDelaun srcVert( 1, 1 + nbP ); + BRepMesh_Vertex v( 0, 0, BRepMesh_Frontier ); + for ( size_t iW = 0; iW < srcWires.size(); ++iW ) + { + const UVPtStructVec& srcPnt = srcWires[iW]->GetUVPtStruct(); + for ( int i = 0, nb = srcPnt.size() - 1; i < nb; ++i, ++iP ) + { + _bndSrcNodes[ iP ] = srcPnt[i].node; + srcPnt[i].node->setIsMarked( true ); + + v.ChangeCoord() = srcPnt[i].UV().Multiplied( _scale ); + srcVert( iP ) = v; + } + } + // triangulate the srcFace in 2D + BRepMesh_Delaun delauney( srcVert ); + _triaDS = delauney.Result(); + } + + //================================================================================ + /*! + * \brief Move non-marked target nodes + * \param [in,out] tgtHelper - helper + * \param [in] tgtWires - boundary nodes of the target FACE; must be in the + * same order as the nodes in srcWires given in the constructor + * \param [in] src2tgtNodes - map of src -> tgt nodes + * \param [in] moveAll - to move all nodes; if \c false, move only non-marked nodes + * \return bool - Ok or not + */ + //================================================================================ + + bool Morph::Perform(SMESH_MesherHelper& tgtHelper, + const TSideVector& tgtWires, + Handle(ShapeAnalysis_Surface) tgtSurface, + const TNodeNodeMap& src2tgtNodes, + const bool moveAll) + { + // get tgt boundary points corresponding to _bndSrcNodes + size_t nbP = 0; + for ( size_t iW = 0; iW < tgtWires.size(); ++iW ) + nbP += tgtWires[iW]->NbPoints() - 1; // 1st and last points coincide + if ( nbP != _bndSrcNodes.size() - 1 ) + return false; + + BRepMesh::Array1OfVertexOfDelaun tgtVert( 1, 1 + nbP ); + BRepMesh_Vertex v( 0, 0, BRepMesh_Frontier ); + for ( size_t iW = 0, iP = 1; iW < tgtWires.size(); ++iW ) + { + const UVPtStructVec& tgtPnt = tgtWires[iW]->GetUVPtStruct(); + for ( int i = 0, nb = tgtPnt.size() - 1; i < nb; ++i, ++iP ) + { + v.ChangeCoord() = tgtPnt[i].UV().Multiplied( _scale ); + tgtVert( iP ) = v; + } + } + + const TopoDS_Face& srcFace = TopoDS::Face( _srcSubMesh->GetSubShape() ); + const int srcFaceID = _srcSubMesh->GetId(); + SMESHDS_Mesh* tgtMesh = tgtHelper.GetMeshDS(); + const SMDS_MeshNode *srcNode, *tgtNode; + const BRepMesh_Triangle *bmTria; + + // initialize a queue of nodes with starting triangles + TNodeTriaList noTriQueue; + size_t iBndSrcN = 1; + for ( ; iBndSrcN < _bndSrcNodes.size() && noTriQueue.empty(); ++iBndSrcN ) + { + // get a triangle + const BRepMesh::ListOfInteger & linkIds = _triaDS->LinksConnectedTo( iBndSrcN ); + const BRepMesh_PairOfIndex & triaIds = _triaDS->ElementsConnectedTo( linkIds.First() ); + const BRepMesh_Triangle& tria = _triaDS->GetElement( triaIds.Index(1) ); + + addCloseNodes( _bndSrcNodes[ iBndSrcN ], &tria, srcFaceID, noTriQueue ); + } + + // un-mark internal src nodes; later we will mark moved nodes + int nbSrcNodes = 0; + SMDS_NodeIteratorPtr nIt = _srcSubMesh->GetSubMeshDS()->GetNodes(); + if ( !nIt || !nIt->more() ) return true; + if ( moveAll ) + { + nbSrcNodes = _srcSubMesh->GetSubMeshDS()->NbNodes(); + while ( nIt->more() ) + nIt->next()->setIsMarked( false ); + } + else + { + while ( nIt->more() ) + nbSrcNodes += int( !nIt->next()->isMarked() ); + } + + // Move tgt nodes + + double bc[3]; // barycentric coordinates + int nodeIDs[3]; + bool checkUV = true; + const SMDS_FacePosition* pos; + + while ( nbSrcNodes > 0 ) + { + while ( !noTriQueue.empty() ) + { + srcNode = noTriQueue.front().first; + bmTria = noTriQueue.front().second; + noTriQueue.pop_front(); + if ( srcNode->isMarked() ) + continue; + --nbSrcNodes; + srcNode->setIsMarked( true ); + + // find a delauney triangle containing the src node + gp_XY uv = tgtHelper.GetNodeUV( srcFace, srcNode, NULL, &checkUV ); + uv *= _scale; + bmTria = findTriangle( uv, bmTria, _triaDS, bc ); + if ( !bmTria ) + continue; + + // compute new coordinates for a corresponding tgt node + gp_XY uvNew( 0., 0. ), nodeUV; + _triaDS->ElementNodes( *bmTria, nodeIDs ); + for ( int i = 0; i < 3; ++i ) + uvNew += bc[i] * tgtVert( nodeIDs[i]).Coord(); + uvNew.SetCoord( uvNew.X() / _scale.X(), uvNew.Y() / _scale.Y() ); + gp_Pnt xyz = tgtSurface->Value( uvNew ); + + // find and move tgt node + TNodeNodeMap::const_iterator n2n = src2tgtNodes.find( srcNode ); + if ( n2n == src2tgtNodes.end() ) continue; + tgtNode = n2n->second; + tgtMesh->MoveNode( tgtNode, xyz.X(), xyz.Y(), xyz.Z() ); + + if (( pos = dynamic_cast< const SMDS_FacePosition* >( tgtNode->GetPosition() ))) + const_cast( pos )->SetParameters( uvNew.X(), uvNew.Y() ); + + addCloseNodes( srcNode, bmTria, srcFaceID, noTriQueue ); + } + + if ( nbSrcNodes > 0 ) + { + // assure that all src nodes are visited + for ( ; iBndSrcN < _bndSrcNodes.size() && noTriQueue.empty(); ++iBndSrcN ) + { + const BRepMesh::ListOfInteger & linkIds = _triaDS->LinksConnectedTo( iBndSrcN ); + const BRepMesh_PairOfIndex & triaIds = _triaDS->ElementsConnectedTo( linkIds.First() ); + const BRepMesh_Triangle& tria = _triaDS->GetElement( triaIds.Index(1) ); + addCloseNodes( _bndSrcNodes[ iBndSrcN ], &tria, srcFaceID, noTriQueue ); + } + if ( noTriQueue.empty() ) + { + SMDS_NodeIteratorPtr nIt = _srcSubMesh->GetSubMeshDS()->GetNodes(); + while ( nIt->more() ) + { + srcNode = nIt->next(); + if ( !srcNode->isMarked() ) + noTriQueue.push_back( make_pair( srcNode, bmTria )); + } + } + } + } + + return true; + + } // Morph::Perform + +gp_XY Morph::GetBndUV(const int iNode) const +{ + return _triaDS->GetNode( iNode ).Coord(); + } + + +} // namespace StdMeshers_ProjectionUtils diff --git a/src/StdMeshers/StdMeshers_ProjectionUtils.hxx b/src/StdMeshers/StdMeshers_ProjectionUtils.hxx index ced8febbd..d44e477f0 100644 --- a/src/StdMeshers/StdMeshers_ProjectionUtils.hxx +++ b/src/StdMeshers/StdMeshers_ProjectionUtils.hxx @@ -30,11 +30,14 @@ #include "SMESH_StdMeshers.hxx" +#include "StdMeshers_FaceSide.hxx" #include "SMDS_MeshElement.hxx" +#include +#include #include -#include #include +#include #include #include #include @@ -136,6 +139,33 @@ namespace StdMeshers_ProjectionUtils bool Invert(); }; + /*! + * \brief Morph mesh on the target FACE to lie within FACE boundary w/o distortion + */ + class Morph + { + std::vector< const SMDS_MeshNode* > _bndSrcNodes; + Handle(BRepMesh_DataStructureOfDelaun) _triaDS; + SMESH_subMesh* _srcSubMesh; + bool _moveAll; + gp_XY _scale; + public: + + Morph(const TSideVector& srcWires); + + bool Perform(SMESH_MesherHelper& tgtHelper, + const TSideVector& tgtWires, + Handle(ShapeAnalysis_Surface) tgtSurface, + const TNodeNodeMap& src2tgtNodes, + const bool moveAll); + + // return source boundary nodes. 0-th node is zero + const std::vector< const SMDS_MeshNode* >& GetBndNodes() const { return _bndSrcNodes; } + + // return UV of the i-th source boundary node + gp_XY GetBndUV(const int iNode) const; + }; + /*! * \brief Looks for association of all sub-shapes of two shapes * \param theShape1 - shape 1 diff --git a/src/StdMeshers/StdMeshers_Projection_2D.cxx b/src/StdMeshers/StdMeshers_Projection_2D.cxx index eeb20f009..4063153dc 100644 --- a/src/StdMeshers/StdMeshers_Projection_2D.cxx +++ b/src/StdMeshers/StdMeshers_Projection_2D.cxx @@ -1179,240 +1179,6 @@ namespace { return true; } - typedef list< pair< const SMDS_MeshNode*, const BRepMesh_Triangle* > > TNodeTriaList; - - //================================================================================ - /*! - * \brief Add in-FACE nodes surrounding a given node to a queue - */ - //================================================================================ - - void addCloseNodes( const SMDS_MeshNode* srcNode, - const BRepMesh_Triangle* bmTria, - const int srcFaceID, - TNodeTriaList & noTriQueue ) - { - // find in-FACE nodes - SMDS_ElemIteratorPtr elems = srcNode->GetInverseElementIterator(SMDSAbs_Face); - while ( elems->more() ) - { - const SMDS_MeshElement* elem = elems->next(); - if ( elem->getshapeId() == srcFaceID ) - { - for ( int i = 0, nb = elem->NbNodes(); i < nb; ++i ) - { - const SMDS_MeshNode* n = elem->GetNode( i ); - if ( !n->isMarked() ) - noTriQueue.push_back( make_pair( n, bmTria )); - } - } - } - } - - //================================================================================ - /*! - * \brief Find a delauney triangle containing a given 2D point and return - * barycentric coordinates within the found triangle - */ - //================================================================================ - - const BRepMesh_Triangle* findTriangle( const gp_XY& uv, - const BRepMesh_Triangle* bmTria, - Handle(BRepMesh_DataStructureOfDelaun)& triaDS, - double bc[3] ) - { - int nodeIDs[3]; - gp_XY nodeUVs[3]; - int linkIDs[3]; - Standard_Boolean ori[3]; - - while ( bmTria ) - { - // check bmTria - - triaDS->ElementNodes( *bmTria, nodeIDs ); - nodeUVs[0] = triaDS->GetNode( nodeIDs[0] ).Coord(); - nodeUVs[1] = triaDS->GetNode( nodeIDs[1] ).Coord(); - nodeUVs[2] = triaDS->GetNode( nodeIDs[2] ).Coord(); - - SMESH_MeshAlgos::GetBarycentricCoords( uv, - nodeUVs[0], nodeUVs[1], nodeUVs[2], - bc[0], bc[1] ); - if ( bc[0] >= 0 && bc[1] >= 0 && bc[0] + bc[1] <= 1 ) - { - bc[2] = 1 - bc[0] - bc[1]; - return bmTria; - } - - // look for a neighbor triangle, which is adjacent to a link intersected - // by a segment( triangle center -> uv ) - - gp_XY gc = ( nodeUVs[0] + nodeUVs[1] + nodeUVs[2] ) / 3.; - gp_XY seg = uv - gc; - - bmTria->Edges( linkIDs, ori ); - int triaID = triaDS->IndexOf( *bmTria ); - bmTria = 0; - - for ( int i = 0; i < 3; ++i ) - { - const BRepMesh_PairOfIndex & triIDs = triaDS->ElementsConnectedTo( linkIDs[i] ); - if ( triIDs.Extent() < 2 ) - continue; // no neighbor triangle - - // check if a link intersects gc2uv - const BRepMesh_Edge & link = triaDS->GetLink( linkIDs[i] ); - const BRepMesh_Vertex & n1 = triaDS->GetNode( link.FirstNode() ); - const BRepMesh_Vertex & n2 = triaDS->GetNode( link.LastNode() ); - gp_XY uv1 = n1.Coord(); - gp_XY lin = n2.Coord() - uv1; // link direction - - double crossSegLin = seg ^ lin; - if ( Abs( crossSegLin ) < std::numeric_limits::min() ) - continue; // parallel - - double uSeg = ( uv1 - gc ) ^ lin / crossSegLin; - if ( 0. <= uSeg && uSeg <= 1. ) - { - bmTria = & triaDS->GetElement( triIDs.Index( 1 + ( triIDs.Index(1) == triaID ))); - break; - } - } - } - return bmTria; - } - - //================================================================================ - /*! - * \brief Morph mesh on the target face to lie within FACE boundary w/o distortion - * - * algo: - * - make a CDT on the src FACE - * - find a triangle containing a src node and get its barycentric coordinates - * - move the node to a point with the same barycentric coordinates in a corresponding - * tgt triangle - */ - //================================================================================ - - bool morph( SMESH_MesherHelper& tgtHelper, - const TopoDS_Face& tgtFace, - const TopoDS_Face& srcFace, - const TSideVector& tgtWires, - const TSideVector& srcWires, - const TAssocTool::TNodeNodeMap& src2tgtNodes ) - { - if ( srcWires.size() != tgtWires.size() ) return false; - if ( srcWires.size() == 1 ) return false; // tmp - - // count boundary points - int iP = 1, nbP = 0; - for ( size_t iW = 0; iW < srcWires.size(); ++iW ) - nbP += srcWires[iW]->NbPoints() - 1; // 1st and last points coincide - - // fill boundary points - BRepMesh::Array1OfVertexOfDelaun srcVert( 1, 1 + nbP ), tgtVert( 1, 1 + nbP ); - vector< const SMDS_MeshNode* > bndSrcNodes( nbP + 1 ); bndSrcNodes[0] = 0; - BRepMesh_Vertex v( 0, 0, BRepMesh_Frontier ); - for ( size_t iW = 0; iW < srcWires.size(); ++iW ) - { - const UVPtStructVec& srcPnt = srcWires[iW]->GetUVPtStruct(); - const UVPtStructVec& tgtPnt = tgtWires[iW]->GetUVPtStruct(); - if ( srcPnt.size() != tgtPnt.size() ) return false; - - for ( int i = 0, nb = srcPnt.size() - 1; i < nb; ++i, ++iP ) - { - bndSrcNodes[ iP ] = srcPnt[i].node; - srcPnt[i].node->setIsMarked( true ); - - v.ChangeCoord() = srcPnt[i].UV(); - srcVert( iP ) = v; - v.ChangeCoord() = tgtPnt[i].UV(); - tgtVert( iP ) = v; - } - } - // triangulate the srcFace in 2D - BRepMesh_Delaun delauney( srcVert ); - Handle(BRepMesh_DataStructureOfDelaun) triaDS = delauney.Result(); - - Handle(ShapeAnalysis_Surface) tgtSurface = tgtHelper.GetSurface( tgtFace ); - SMESHDS_Mesh* srcMesh = srcWires[0]->GetMesh()->GetMeshDS(); - SMESHDS_Mesh* tgtMesh = tgtHelper.GetMeshDS(); - const SMDS_MeshNode *srcNode, *tgtNode; - const BRepMesh_Triangle *bmTria; - - // un-mark internal src nodes; later we will mark moved nodes - SMDS_NodeIteratorPtr nIt = srcMesh->MeshElements( srcFace )->GetNodes(); - if ( !nIt || !nIt->more() ) return true; - while ( nIt->more() ) - ( srcNode = nIt->next() )->setIsMarked( false ); - - // initialize a queue of nodes with starting triangles - const int srcFaceID = srcNode->getshapeId(); - TNodeTriaList noTriQueue; - size_t iBndSrcN = 1; - for ( ; iBndSrcN < bndSrcNodes.size() && noTriQueue.empty(); ++iBndSrcN ) - { - // get a triangle - const BRepMesh::ListOfInteger & linkIds = triaDS->LinksConnectedTo( iBndSrcN ); - const BRepMesh_PairOfIndex & triaIds = triaDS->ElementsConnectedTo( linkIds.First() ); - const BRepMesh_Triangle& tria = triaDS->GetElement( triaIds.Index(1) ); - - addCloseNodes( bndSrcNodes[ iBndSrcN ], &tria, srcFaceID, noTriQueue ); - } - - // Move tgt nodes - - double bc[3]; // barycentric coordinates - int nodeIDs[3]; - bool checkUV = true; - const SMDS_FacePosition* pos; - - while ( !noTriQueue.empty() ) - { - srcNode = noTriQueue.front().first; - bmTria = noTriQueue.front().second; - noTriQueue.pop_front(); - if ( srcNode->isMarked() ) - continue; - srcNode->setIsMarked( true ); - - // find a delauney triangle containing the src node - gp_XY uv = tgtHelper.GetNodeUV( srcFace, srcNode, NULL, &checkUV ); - bmTria = findTriangle( uv, bmTria, triaDS, bc ); - if ( !bmTria ) - continue; - - // compute new coordinates for a corresponding tgt node - gp_XY uvNew( 0., 0. ), nodeUV; - triaDS->ElementNodes( *bmTria, nodeIDs ); - for ( int i = 0; i < 3; ++i ) - uvNew += bc[i] * tgtVert( nodeIDs[i]).Coord(); - gp_Pnt xyz = tgtSurface->Value( uvNew ); - - // find and move tgt node - TAssocTool::TNodeNodeMap::const_iterator n2n = src2tgtNodes.find( srcNode ); - if ( n2n == src2tgtNodes.end() ) continue; - tgtNode = n2n->second; - tgtMesh->MoveNode( tgtNode, xyz.X(), xyz.Y(), xyz.Z() ); - - if (( pos = dynamic_cast< const SMDS_FacePosition* >( tgtNode->GetPosition() ))) - const_cast( pos )->SetParameters( uvNew.X(), uvNew.Y() ); - - addCloseNodes( srcNode, bmTria, srcFaceID, noTriQueue ); - - // assure that all src nodes are visited - for ( ; iBndSrcN < bndSrcNodes.size() && noTriQueue.empty(); ++iBndSrcN ) - { - const BRepMesh::ListOfInteger & linkIds = triaDS->LinksConnectedTo( iBndSrcN ); - const BRepMesh_PairOfIndex & triaIds = triaDS->ElementsConnectedTo( linkIds.First() ); - const BRepMesh_Triangle& tria = triaDS->GetElement( triaIds.Index(1) ); - addCloseNodes( bndSrcNodes[ iBndSrcN ], &tria, srcFaceID, noTriQueue ); - } - } - - return true; - } - //======================================================================= /* * Set initial association of VERTEXes for the case of projection @@ -1933,7 +1699,9 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& // ---------------------------------------------------------------- if ( helper.IsDistorted2D( tgtSubMesh, /*checkUV=*/false, &helper )) { - morph( helper, tgtFace, srcFace, tgtWires, srcWires, _src2tgtNodes ); + TAssocTool::Morph morph( srcWires ); + morph.Perform( helper, tgtWires, helper.GetSurface( tgtFace ), + _src2tgtNodes, /*moveAll=*/true ); if ( !fixDistortedFaces( helper, tgtWires )) return error("Invalid mesh generated"); diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index a6530b924..8212b0a8b 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -95,7 +95,7 @@ #include #ifdef _DEBUG_ -//#define __myDEBUG +#define __myDEBUG //#define __NOT_INVALIDATE_BAD_SMOOTH //#define __NODES_AT_POS #endif @@ -427,20 +427,21 @@ namespace VISCOUS_3D enum EFlags { TO_SMOOTH = 0x0000001, MOVED = 0x0000002, // set by _neibors[i]->SetNewLength() - SMOOTHED = 0x0000004, // set by this->Smooth() + SMOOTHED = 0x0000004, // set by _LayerEdge::Smooth() DIFFICULT = 0x0000008, // near concave VERTEX ON_CONCAVE_FACE = 0x0000010, BLOCKED = 0x0000020, // not to inflate any more INTERSECTED = 0x0000040, // close intersection with a face found NORMAL_UPDATED = 0x0000080, - MARKED = 0x0000100, // local usage - MULTI_NORMAL = 0x0000200, // a normal is invisible by some of surrounding faces - NEAR_BOUNDARY = 0x0000400, // is near FACE boundary forcing smooth - SMOOTHED_C1 = 0x0000800, // is on _eosC1 - DISTORTED = 0x0001000, // was bad before smoothing - RISKY_SWOL = 0x0002000, // SWOL is parallel to a source FACE - SHRUNK = 0x0004000, // target node reached a tgt position while shrink() - UNUSED_FLAG = 0x0100000 // to add use flags after + UPD_NORMAL_CONV = 0x0000100, // to update normal on boundary of concave FACE + MARKED = 0x0000200, // local usage + MULTI_NORMAL = 0x0000400, // a normal is invisible by some of surrounding faces + NEAR_BOUNDARY = 0x0000800, // is near FACE boundary forcing smooth + SMOOTHED_C1 = 0x0001000, // is on _eosC1 + DISTORTED = 0x0002000, // was bad before smoothing + RISKY_SWOL = 0x0004000, // SWOL is parallel to a source FACE + SHRUNK = 0x0008000, // target node reached a tgt position while shrink() + UNUSED_FLAG = 0x0100000 // to add user flags after }; bool Is ( int flag ) const { return _flags & flag; } void Set ( int flag ) { _flags |= flag; } @@ -497,7 +498,7 @@ namespace VISCOUS_3D const gp_XYZ& PrevPos() const { return _pos[ _pos.size() - 2 ]; } gp_XYZ PrevCheckPos( _EdgesOnShape* eos=0 ) const; gp_Ax1 LastSegment(double& segLen, _EdgesOnShape& eos) const; - gp_XY LastUV( const TopoDS_Face& F, _EdgesOnShape& eos ) const; + gp_XY LastUV( const TopoDS_Face& F, _EdgesOnShape& eos, int which=-1 ) const; bool IsOnEdge() const { return _2neibors; } gp_XYZ Copy( _LayerEdge& other, _EdgesOnShape& eos, SMESH_MesherHelper& helper ); void SetCosin( double cosin ); @@ -664,6 +665,8 @@ namespace VISCOUS_3D _SolidData* _data; // parent SOLID + _LayerEdge* operator[](size_t i) const { return (_LayerEdge*) _edges[i]; } + size_t size() const { return _edges.size(); } TopAbs_ShapeEnum ShapeType() const { return _shape.IsNull() ? TopAbs_SHAPE : _shape.ShapeType(); } TopAbs_ShapeEnum SWOLType() const @@ -678,7 +681,7 @@ namespace VISCOUS_3D //-------------------------------------------------------------------------------- /*! - * \brief Convex FACE whose radius of curvature is less than the thickness of + * \brief Convex FACE whose radius of curvature is less than the thickness of * layers. It is used to detect distortion of prisms based on a convex * FACE and to update normals to enable further increasing the thickness */ @@ -692,7 +695,14 @@ namespace VISCOUS_3D // map a sub-shape to _SolidData::_edgesOnShape map< TGeomID, _EdgesOnShape* > _subIdToEOS; + bool _isTooCurved; bool _normalsFixed; + bool _normalsFixedOnBorders; // used in putOnOffsetSurface() + + double GetMaxCurvature( _SolidData& data, + _EdgesOnShape& eof, + BRepLProp_SLProps& surfProp, + SMESH_MesherHelper& helper); bool GetCenterOfCurvature( _LayerEdge* ledge, BRepLProp_SLProps& surfProp, @@ -757,8 +767,6 @@ namespace VISCOUS_3D int _nbShapesToSmooth; - //map< TGeomID,Handle(Geom_Curve)> _edge2curve; - vector< _CollisionEdges > _collisionEdges; set< TGeomID > _concaveFaces; @@ -870,6 +878,7 @@ namespace VISCOUS_3D const gp_XY& uvToFix, const double refSign ); }; + struct PyDump; //-------------------------------------------------------------------------------- /*! * \brief Builder of viscous layers @@ -942,8 +951,11 @@ namespace VISCOUS_3D void makeOffsetSurface( _EdgesOnShape& eos, SMESH_MesherHelper& ); void putOnOffsetSurface( _EdgesOnShape& eos, int infStep, vector< _EdgesOnShape* >& eosC1, - int smooStep=0, bool moveAll=false ); + int smooStep=0, int moveAll=false ); void findCollisionEdges( _SolidData& data, SMESH_MesherHelper& helper ); + void findEdgesToUpdateNormalNearConvexFace( _ConvexFace & convFace, + _SolidData& data, + SMESH_MesherHelper& helper ); void limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelper& helper ); void limitMaxLenByCurvature( _LayerEdge* e1, _LayerEdge* e2, _EdgesOnShape& eos1, _EdgesOnShape& eos2, @@ -987,6 +999,7 @@ namespace VISCOUS_3D TopTools_MapOfShape _shrinkedFaces; int _tmpFaceID; + PyDump* _pyDump; }; //-------------------------------------------------------------------------------- /*! @@ -1034,6 +1047,7 @@ namespace VISCOUS_3D size_t _iSeg[2]; // index of segment where extreme tgt node is projected _EdgesOnShape& _eos; double _curveLen; // length of the EDGE + std::pair _eToSmooth[2]; // indices of _LayerEdge's in _eos static Handle(Geom_Curve) CurveForSmooth( const TopoDS_Edge& E, _EdgesOnShape& eos, @@ -1047,31 +1061,24 @@ namespace VISCOUS_3D bool Perform(_SolidData& data, Handle(ShapeAnalysis_Surface)& surface, const TopoDS_Face& F, - SMESH_MesherHelper& helper ) - { - if ( _leParams.empty() || ( !isAnalytic() && _offPoints.empty() )) - prepare( data ); + SMESH_MesherHelper& helper ); - if ( isAnalytic() ) - return smoothAnalyticEdge( data, surface, F, helper ); - else - return smoothComplexEdge ( data, surface, F, helper ); - } void prepare(_SolidData& data ); + void findEdgesToSmooth(); + + bool isToSmooth( int iE ); + bool smoothAnalyticEdge( _SolidData& data, Handle(ShapeAnalysis_Surface)& surface, const TopoDS_Face& F, SMESH_MesherHelper& helper); - bool smoothComplexEdge( _SolidData& data, Handle(ShapeAnalysis_Surface)& surface, const TopoDS_Face& F, SMESH_MesherHelper& helper); - gp_XYZ getNormalNormal( const gp_XYZ & normal, const gp_XYZ& edgeDir); - _LayerEdge* getLEdgeOnV( bool is2nd ) { return _eos._edges[ is2nd ? _eos._edges.size()-1 : 0 ]->_2neibors->_edges[ is2nd ]; @@ -1663,14 +1670,15 @@ namespace VISCOUS_3D // HOWTO use: run python commands written in a console to see // construction steps of viscous layers #ifdef __myDEBUG - ofstream* py; - int theNbPyFunc; - struct PyDump { + ostream* py; + int theNbPyFunc; + struct PyDump + { PyDump(SMESH_Mesh& m) { int tag = 3 + m.GetId(); const char* fname = "/tmp/viscous.py"; cout << "execfile('"< ostream & operator<<( const T &anything ) { return *this ; } + }; + void Pause() { py = &_mystream; } + void Resume() { py = _pyStream; } + MyStream _mystream; + ostream* _pyStream; }; #define dumpFunction(f) { _dumpFunction(f, __LINE__);} #define dumpMove(n) { _dumpMove(n, __LINE__);} @@ -1710,7 +1726,7 @@ namespace VISCOUS_3D #else - struct PyDump { PyDump(SMESH_Mesh&) {} void Finish() {} }; + struct PyDump { PyDump(SMESH_Mesh&) {} void Finish() {} void Pause() {} void Resume() {} }; #define dumpFunction(f) f #define dumpMove(n) #define dumpMoveComm(n,txt) @@ -1852,6 +1868,7 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh& theMesh, return SMESH_ComputeErrorPtr(); // everything already computed PyDump debugDump( theMesh ); + _pyDump = &debugDump; // TODO: ignore already computed SOLIDs if ( !findSolidsWithLayers()) @@ -2777,8 +2794,6 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) { SMESH_MesherHelper helper( *_mesh ); - const int nbTestPnt = 5; // on a FACE sub-shape - BRepLProp_SLProps surfProp( 2, 1e-6 ); data._convexFaces.clear(); @@ -2790,58 +2805,27 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) continue; TopoDS_Face F = TopoDS::Face( eof._shape ); - SMESH_subMesh * sm = eof._subMesh; const TGeomID faceID = eof._shapeID; BRepAdaptor_Surface surface( F, false ); surfProp.SetSurface( surface ); - bool isTooCurved = false; - _ConvexFace cnvFace; - const double oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. ); - SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true); - while ( smIt->more() ) - { - sm = smIt->next(); - const TGeomID subID = sm->GetId(); - // find _LayerEdge's of a sub-shape - _EdgesOnShape* eos; - if (( eos = data.GetShapeEdges( subID ))) - cnvFace._subIdToEOS.insert( make_pair( subID, eos )); - else - continue; - // check concavity and curvature and limit data._stepSize - const double minCurvature = - 1. / ( eos->_hyp.GetTotalThickness() * ( 1 + theThickToIntersection )); - size_t iStep = Max( 1, eos->_edges.size() / nbTestPnt ); - for ( size_t i = 0; i < eos->_edges.size(); i += iStep ) - { - gp_XY uv = helper.GetNodeUV( F, eos->_edges[ i ]->_nodes[0] ); - surfProp.SetParameters( uv.X(), uv.Y() ); - if ( !surfProp.IsCurvatureDefined() ) - continue; - if ( surfProp.MaxCurvature() * oriFactor > minCurvature ) - { - limitStepSize( data, 0.9 / surfProp.MaxCurvature() * oriFactor ); - isTooCurved = true; - } - if ( surfProp.MinCurvature() * oriFactor > minCurvature ) - { - limitStepSize( data, 0.9 / surfProp.MinCurvature() * oriFactor ); - isTooCurved = true; - } - } - } // loop on sub-shapes of the FACE + cnvFace._face = F; + cnvFace._normalsFixed = false; + cnvFace._isTooCurved = false; - if ( !isTooCurved ) continue; + double maxCurvature = cnvFace.GetMaxCurvature( data, eof, surfProp, helper ); + if ( maxCurvature > 0 ) + { + limitStepSize( data, 0.9 / maxCurvature ); + findEdgesToUpdateNormalNearConvexFace( cnvFace, data, helper ); + } + if ( !cnvFace._isTooCurved ) continue; _ConvexFace & convFace = data._convexFaces.insert( make_pair( faceID, cnvFace )).first->second; - convFace._face = F; - convFace._normalsFixed = false; - // skip a closed surface (data._convexFaces is useful anyway) bool isClosedF = false; helper.SetSubShape( F ); @@ -2854,6 +2838,7 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) if ( isClosedF ) { // limit _LayerEdge::_maxLen on the FACE + const double oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. ); const double minCurvature = 1. / ( eof._hyp.GetTotalThickness() * ( 1 + theThickToIntersection )); map< TGeomID, _EdgesOnShape* >::iterator id2eos = cnvFace._subIdToEOS.find( faceID ); @@ -2865,14 +2850,13 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) _LayerEdge* ledge = eos._edges[ i ]; gp_XY uv = helper.GetNodeUV( F, ledge->_nodes[0] ); surfProp.SetParameters( uv.X(), uv.Y() ); - if ( !surfProp.IsCurvatureDefined() ) - continue; - - if ( surfProp.MaxCurvature() * oriFactor > minCurvature ) - ledge->_maxLen = Min( ledge->_maxLen, 1. / surfProp.MaxCurvature() * oriFactor ); - - if ( surfProp.MinCurvature() * oriFactor > minCurvature ) - ledge->_maxLen = Min( ledge->_maxLen, 1. / surfProp.MinCurvature() * oriFactor ); + if ( surfProp.IsCurvatureDefined() ) + { + double curvature = Max( surfProp.MaxCurvature() * oriFactor, + surfProp.MinCurvature() * oriFactor ); + if ( curvature > minCurvature ) + ledge->_maxLen = Min( ledge->_maxLen, 1. / curvature ); + } } } continue; @@ -3049,8 +3033,8 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data ) { eos._edgeSmoother = new _Smoother1D( curve, eos ); - for ( size_t i = 0; i < eos._edges.size(); ++i ) - eos._edges[i]->Set( _LayerEdge::TO_SMOOTH ); + // for ( size_t i = 0; i < eos._edges.size(); ++i ) + // eos._edges[i]->Set( _LayerEdge::TO_SMOOTH ); } } } @@ -3429,11 +3413,12 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, } // find _normal + bool fromVonF = false; if ( useGeometry ) { - bool fromVonF = ( eos.ShapeType() == TopAbs_VERTEX && - eos.SWOLType() == TopAbs_FACE && - totalNbFaces > 1 ); + fromVonF = ( eos.ShapeType() == TopAbs_VERTEX && + eos.SWOLType() == TopAbs_FACE && + totalNbFaces > 1 ); if ( onShrinkShape && !fromVonF ) // one of faces the node is on has no layers { @@ -3535,14 +3520,19 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, break; } case TopAbs_VERTEX: { - //if ( eos.SWOLType() != TopAbs_FACE ) // else _cosin is set by getFaceDir() + if ( fromVonF ) + { + getFaceDir( TopoDS::Face( eos._sWOL ), TopoDS::Vertex( eos._shape ), + node, helper, normOK, &edge._cosin ); + } + else if ( eos.SWOLType() != TopAbs_FACE ) // else _cosin is set by getFaceDir() { TopoDS_Vertex V = TopoDS::Vertex( eos._shape ); gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK ); double angle = inFaceDir.Angle( edge._normal ); // [0,PI] edge._cosin = Cos( angle ); if ( totalNbFaces > 2 || helper.IsSeamShape( node->getshapeId() )) - for ( int iF = totalNbFaces-2; iF >=0; --iF ) + for ( int iF = 1; iF < totalNbFaces; ++iF ) { F = face2Norm[ iF ].first; inFaceDir = getFaceDir( F, V, node, helper, normOK=true ); @@ -4428,11 +4418,14 @@ bool _ViscousBuilder::inflate(_SolidData& data) data._epsilon = data._stepSize * 1e-7; debugMsg( "-- geomSize = " << data._geomSize << ", stepSize = " << data._stepSize ); + _pyDump->Pause(); findCollisionEdges( data, helper ); limitMaxLenByCurvature( data, helper ); + _pyDump->Resume(); + // limit length of _LayerEdge's around MULTI_NORMAL _LayerEdge's for ( size_t i = 0; i < data._edgesOnShape.size(); ++i ) if ( data._edgesOnShape[i].ShapeType() == TopAbs_VERTEX && @@ -4959,7 +4952,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, // ignore intersection of a _LayerEdge based on a _ConvexFace with a face // lying on this _ConvexFace if ( _ConvexFace* convFace = data.GetConvexFace( intFace->getshapeId() )) - if ( convFace->_subIdToEOS.count ( eos._shapeID )) + if ( convFace->_isTooCurved && convFace->_subIdToEOS.count ( eos._shapeID )) continue; // ignore intersection of a _LayerEdge based on a FACE with an element on this FACE @@ -4971,15 +4964,16 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, if ( dist > 0 ) { bool toIgnore = false; - if ( eos._edges[i]->Is( _LayerEdge::TO_SMOOTH )) + if ( eos._toSmooth ) { const TopoDS_Shape& S = getMeshDS()->IndexToShape( intFace->getshapeId() ); if ( !S.IsNull() && S.ShapeType() == TopAbs_FACE ) { - TopExp_Explorer edge( eos._shape, TopAbs_EDGE ); - for ( ; !toIgnore && edge.More(); edge.Next() ) - // is adjacent - has a common EDGE - toIgnore = ( helper.IsSubShape( edge.Current(), S )); + TopExp_Explorer sub( eos._shape, + eos.ShapeType() == TopAbs_FACE ? TopAbs_EDGE : TopAbs_VERTEX ); + for ( ; !toIgnore && sub.More(); sub.Next() ) + // is adjacent - has a common EDGE or VERTEX + toIgnore = ( helper.IsSubShape( sub.Current(), S )); if ( toIgnore ) // check angle between normals { @@ -5264,7 +5258,7 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, int infStep, vector< _EdgesOnShape* >& eosC1, int smooStep, - bool moveAll ) + int moveAll ) { _EdgesOnShape * eof = & eos; if ( eos.ShapeType() != TopAbs_FACE ) // eos is a boundary of C1 FACE, look for the FACE eos @@ -5295,8 +5289,13 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, edge->Unset( _LayerEdge::MARKED ); if ( edge->Is( _LayerEdge::BLOCKED ) || !edge->_curvature ) continue; - if ( !moveAll && !edge->Is( _LayerEdge::MOVED )) + if ( moveAll == _LayerEdge::UPD_NORMAL_CONV ) + { + if ( !edge->Is( _LayerEdge::UPD_NORMAL_CONV )) continue; + } + else if ( !moveAll && !edge->Is( _LayerEdge::MOVED )) + continue; int nbBlockedAround = 0; for ( size_t iN = 0; iN < edge->_neibors.size(); ++iN ) @@ -5306,7 +5305,7 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, gp_Pnt tgtP = SMESH_TNodeXYZ( edge->_nodes.back() ); gp_Pnt2d uv = eof->_offsetSurf->NextValueOfUV( edge->_curvature->_uv, tgtP, preci ); - if ( eof->_offsetSurf->Gap() > edge->_len ) continue; // NextValueOfUV() bug + if ( eof->_offsetSurf->Gap() > edge->_len ) continue; // NextValueOfUV() bug edge->_curvature->_uv = uv; if ( eof->_offsetSurf->Gap() < 10 * preci ) continue; // same pos @@ -5325,9 +5324,15 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, edge->_pos.back() = newP; edge->Set( _LayerEdge::MARKED ); + if ( moveAll == _LayerEdge::UPD_NORMAL_CONV ) + { + edge->_normal = ( newP - prevP ).Normalized(); + } } } + + #ifdef _DEBUG_ // dumpMove() for debug size_t i = 0; @@ -5336,7 +5341,7 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, break; if ( i < eos._edges.size() ) { - dumpFunction(SMESH_Comment("putOnOffsetSurface_F") << eos._shapeID + dumpFunction(SMESH_Comment("putOnOffsetSurface_S") << eos._shapeID << "_InfStep" << infStep << "_" << smooStep ); for ( ; i < eos._edges.size(); ++i ) { @@ -5346,6 +5351,26 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, dumpFunctionEnd(); } #endif + + _ConvexFace* cnvFace; + if ( moveAll != _LayerEdge::UPD_NORMAL_CONV && + eos.ShapeType() == TopAbs_FACE && + (cnvFace = eos.GetData().GetConvexFace( eos._shapeID )) && + !cnvFace->_normalsFixedOnBorders ) + { + // put on the surface nodes built on FACE boundaries + SMESH_subMeshIteratorPtr smIt = eos._subMesh->getDependsOnIterator(/*includeSelf=*/false); + while ( smIt->more() ) + { + SMESH_subMesh* sm = smIt->next(); + _EdgesOnShape* subEOS = eos.GetData().GetShapeEdges( sm->GetId() ); + if ( !subEOS->_sWOL.IsNull() ) continue; + if ( std::find( eosC1.begin(), eosC1.end(), subEOS ) != eosC1.end() ) continue; + + putOnOffsetSurface( *subEOS, infStep, eosC1, smooStep, _LayerEdge::UPD_NORMAL_CONV ); + } + cnvFace->_normalsFixedOnBorders = true; + } } //================================================================================ @@ -5443,6 +5468,106 @@ Handle(Geom_Curve) _Smoother1D::CurveForSmooth( const TopoDS_Edge& E, return Handle(Geom_Curve)(); } +//================================================================================ +/*! + * \brief Smooth edges on EDGE + */ +//================================================================================ + +bool _Smoother1D::Perform(_SolidData& data, + Handle(ShapeAnalysis_Surface)& surface, + const TopoDS_Face& F, + SMESH_MesherHelper& helper ) +{ + if ( _leParams.empty() || ( !isAnalytic() && _offPoints.empty() )) + prepare( data ); + + findEdgesToSmooth(); + if ( isAnalytic() ) + return smoothAnalyticEdge( data, surface, F, helper ); + else + return smoothComplexEdge ( data, surface, F, helper ); +} + +//================================================================================ +/*! + * \brief Find edges to smooth + */ +//================================================================================ + +void _Smoother1D::findEdgesToSmooth() +{ + _LayerEdge* leOnV[2] = { getLEdgeOnV(0), getLEdgeOnV(1) }; + for ( int iEnd = 0; iEnd < 2; ++iEnd ) + if ( leOnV[iEnd]->Is( _LayerEdge::NORMAL_UPDATED )) + _leOnV[iEnd]._cosin = Abs( _edgeDir[iEnd].Normalized() * leOnV[iEnd]->_normal ); + + _eToSmooth[0].first = _eToSmooth[0].second = 0; + + for ( size_t i = 0; i < _eos.size(); ++i ) + { + if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH )) + { + if ( needSmoothing( _leOnV[0]._cosin, _eos[i]->_len, _curveLen * _leParams[i] ) || + isToSmooth( i )) + _eos[i]->Set( _LayerEdge::TO_SMOOTH ); + else + break; + } + _eToSmooth[0].second = i+1; + } + + _eToSmooth[1].first = _eToSmooth[1].second = _eos.size(); + + for ( int i = _eos.size() - 1; i >= _eToSmooth[0].second; --i ) + { + if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH )) + { + if ( needSmoothing( _leOnV[1]._cosin, _eos[i]->_len, _curveLen * ( 1.-_leParams[i] )) || + isToSmooth( i )) + _eos[i]->Set( _LayerEdge::TO_SMOOTH ); + else + break; + } + _eToSmooth[1].first = i; + } +} + +//================================================================================ +/*! + * \brief Check if iE-th _LayerEdge needs smoothing + */ +//================================================================================ + +bool _Smoother1D::isToSmooth( int iE ) +{ + SMESH_NodeXYZ pi( _eos[iE]->_nodes[0] ); + SMESH_NodeXYZ p0( _eos[iE]->_2neibors->srcNode(0) ); + SMESH_NodeXYZ p1( _eos[iE]->_2neibors->srcNode(1) ); + gp_XYZ seg0 = pi - p0; + gp_XYZ seg1 = p1 - pi; + gp_XYZ tangent = seg0 + seg1; + double tangentLen = tangent.Modulus(); + double segMinLen = Min( seg0.Modulus(), seg1.Modulus() ); + if ( tangentLen < std::numeric_limits::min() ) + return false; + tangent /= tangentLen; + + for ( size_t i = 0; i < _eos[iE]->_neibors.size(); ++i ) + { + _LayerEdge* ne = _eos[iE]->_neibors[i]; + if ( !ne->Is( _LayerEdge::TO_SMOOTH ) || + ne->_nodes.size() < 2 || + ne->_nodes[0]->GetPosition()->GetDim() != 2 ) + continue; + gp_XYZ edgeVec = SMESH_NodeXYZ( ne->_nodes.back() ) - SMESH_NodeXYZ( ne->_nodes[0] ); + double proj = edgeVec * tangent; + if ( needSmoothing( 1., proj, segMinLen )) + return true; + } + return false; +} + //================================================================================ /*! * \brief smooth _LayerEdge's on a staight EDGE or circular EDGE @@ -5456,80 +5581,100 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData& data, { if ( !isAnalytic() ) return false; - const size_t iFrom = 0, iTo = _eos._edges.size(); + size_t iFrom = 0, iTo = _eos._edges.size(); if ( _anaCurve->IsKind( STANDARD_TYPE( Geom_Line ))) { if ( F.IsNull() ) // 3D { - SMESH_TNodeXYZ p0 ( _eos._edges[iFrom]->_2neibors->tgtNode(0) ); - SMESH_TNodeXYZ p1 ( _eos._edges[iTo-1]->_2neibors->tgtNode(1) ); SMESH_TNodeXYZ pSrc0( _eos._edges[iFrom]->_2neibors->srcNode(0) ); SMESH_TNodeXYZ pSrc1( _eos._edges[iTo-1]->_2neibors->srcNode(1) ); - gp_XYZ newPos, lineDir = pSrc1 - pSrc0; - _LayerEdge* vLE0 = _eos._edges[iFrom]->_2neibors->_edges[0]; - _LayerEdge* vLE1 = _eos._edges[iTo-1]->_2neibors->_edges[1]; - bool shiftOnly = ( vLE0->Is( _LayerEdge::NORMAL_UPDATED ) || - vLE0->Is( _LayerEdge::BLOCKED ) || - vLE1->Is( _LayerEdge::NORMAL_UPDATED ) || - vLE1->Is( _LayerEdge::BLOCKED )); - for ( size_t i = iFrom; i < iTo; ++i ) - { - _LayerEdge* edge = _eos._edges[i]; - SMDS_MeshNode* tgtNode = const_cast( edge->_nodes.back() ); - newPos = p0 * ( 1. - _leParams[i] ) + p1 * _leParams[i]; - - if ( shiftOnly || edge->Is( _LayerEdge::NORMAL_UPDATED )) - { - gp_XYZ curPos = SMESH_TNodeXYZ ( tgtNode ); - double shift = ( lineDir * ( newPos - pSrc0 ) - - lineDir * ( curPos - pSrc0 )); - newPos = curPos + lineDir * shift / lineDir.SquareModulus(); - } - if ( edge->Is( _LayerEdge::BLOCKED )) + //const gp_XYZ lineDir = pSrc1 - pSrc0; + //_LayerEdge* vLE0 = getLEdgeOnV( 0 ); + //_LayerEdge* vLE1 = getLEdgeOnV( 1 ); + // bool shiftOnly = ( vLE0->Is( _LayerEdge::NORMAL_UPDATED ) || + // vLE0->Is( _LayerEdge::BLOCKED ) || + // vLE1->Is( _LayerEdge::NORMAL_UPDATED ) || + // vLE1->Is( _LayerEdge::BLOCKED )); + for ( int iEnd = 0; iEnd < 2; ++iEnd ) + { + iFrom = _eToSmooth[ iEnd ].first, iTo = _eToSmooth[ iEnd ].second; + if ( iFrom >= iTo ) continue; + SMESH_TNodeXYZ p0( _eos[iFrom]->_2neibors->tgtNode(0) ); + SMESH_TNodeXYZ p1( _eos[iTo-1]->_2neibors->tgtNode(1) ); + double param0 = ( iFrom == 0 ) ? 0. : _leParams[ iFrom-1 ]; + double param1 = _leParams[ iTo ]; + for ( size_t i = iFrom; i < iTo; ++i ) { - SMESH_TNodeXYZ pSrc( edge->_nodes[0] ); - double curThick = pSrc.SquareDistance( tgtNode ); - double newThink = ( pSrc - newPos ).SquareModulus(); - if ( newThink > curThick ) - continue; + _LayerEdge* edge = _eos[i]; + SMDS_MeshNode* tgtNode = const_cast( edge->_nodes.back() ); + double param = ( _leParams[i] - param0 ) / ( param1 - param0 ); + gp_XYZ newPos = p0 * ( 1. - param ) + p1 * param; + + // if ( shiftOnly || edge->Is( _LayerEdge::NORMAL_UPDATED )) + // { + // gp_XYZ curPos = SMESH_TNodeXYZ ( tgtNode ); + // double shift = ( lineDir * ( newPos - pSrc0 ) - + // lineDir * ( curPos - pSrc0 )); + // newPos = curPos + lineDir * shift / lineDir.SquareModulus(); + // } + if ( edge->Is( _LayerEdge::BLOCKED )) + { + SMESH_TNodeXYZ pSrc( edge->_nodes[0] ); + double curThick = pSrc.SquareDistance( tgtNode ); + double newThink = ( pSrc - newPos ).SquareModulus(); + if ( newThink > curThick ) + continue; + } + edge->_pos.back() = newPos; + tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() ); + dumpMove( tgtNode ); } - edge->_pos.back() = newPos; - tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() ); - dumpMove( tgtNode ); } } else // 2D { - _LayerEdge* e0 = getLEdgeOnV( 0 ); - _LayerEdge* e1 = getLEdgeOnV( 1 ); - gp_XY uv0 = e0->LastUV( F, *data.GetShapeEdges( e0 )); - gp_XY uv1 = e1->LastUV( F, *data.GetShapeEdges( e1 )); - if ( e0->_nodes.back() == e1->_nodes.back() ) // closed edge + _LayerEdge* eV0 = getLEdgeOnV( 0 ); + _LayerEdge* eV1 = getLEdgeOnV( 1 ); + gp_XY uvV0 = eV0->LastUV( F, *data.GetShapeEdges( eV0 )); + gp_XY uvV1 = eV1->LastUV( F, *data.GetShapeEdges( eV1 )); + if ( eV0->_nodes.back() == eV1->_nodes.back() ) // closed edge { int iPeriodic = helper.GetPeriodicIndex(); if ( iPeriodic == 1 || iPeriodic == 2 ) { - uv1.SetCoord( iPeriodic, helper.GetOtherParam( uv1.Coord( iPeriodic ))); - if ( uv0.Coord( iPeriodic ) > uv1.Coord( iPeriodic )) - std::swap( uv0, uv1 ); + uvV1.SetCoord( iPeriodic, helper.GetOtherParam( uvV1.Coord( iPeriodic ))); + if ( uvV0.Coord( iPeriodic ) > uvV1.Coord( iPeriodic )) + std::swap( uvV0, uvV1 ); } } - const gp_XY rangeUV = uv1 - uv0; - for ( size_t i = iFrom; i < iTo; ++i ) - { - if ( _eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue; - gp_XY newUV = uv0 + _leParams[i] * rangeUV; - _eos._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 ); - - gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() ); - SMDS_MeshNode* tgtNode = const_cast( _eos._edges[i]->_nodes.back() ); - tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() ); - dumpMove( tgtNode ); + for ( int iEnd = 0; iEnd < 2; ++iEnd ) + { + iFrom = _eToSmooth[ iEnd ].first, iTo = _eToSmooth[ iEnd ].second; + if ( iFrom >= iTo ) continue; + _LayerEdge* e0 = _eos[iFrom]->_2neibors->_edges[0]; + _LayerEdge* e1 = _eos[iTo-1]->_2neibors->_edges[1]; + gp_XY uv0 = ( e0 == eV0 ) ? uvV0 : e0->LastUV( F, _eos ); + gp_XY uv1 = ( e1 == eV1 ) ? uvV1 : e1->LastUV( F, _eos ); + double param0 = ( iFrom == 0 ) ? 0. : _leParams[ iFrom-1 ]; + double param1 = _leParams[ iTo ]; + const gp_XY rangeUV = uv1 - uv0; + for ( size_t i = iFrom; i < iTo; ++i ) + { + if ( _eos[i]->Is( _LayerEdge::BLOCKED )) continue; + double param = ( _leParams[i] - param0 ) / ( param1 - param0 ); + gp_XY newUV = uv0 + param * rangeUV; + _eos[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 ); + + gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() ); + SMDS_MeshNode* tgtNode = const_cast( _eos[i]->_nodes.back() ); + tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() ); + dumpMove( tgtNode ); - SMDS_FacePosition* pos = static_cast( tgtNode->GetPosition() ); - pos->SetUParameter( newUV.X() ); - pos->SetVParameter( newUV.Y() ); + SMDS_FacePosition* pos = static_cast( tgtNode->GetPosition() ); + pos->SetUParameter( newUV.X() ); + pos->SetVParameter( newUV.Y() ); + } } } return true; @@ -5568,9 +5713,10 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData& data, if ( uLast < 0 ) uLast += 2 * M_PI; - for ( size_t i = iFrom; i < iTo; ++i ) + for ( size_t i = 0; i < _eos.size(); ++i ) { - if ( _eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue; + if ( _eos[i]->Is( _LayerEdge::BLOCKED )) continue; + //if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH )) continue; double u = uLast * _leParams[i]; gp_Pnt p = ElCLib::Value( u, newCirc ); _eos._edges[i]->_pos.back() = p.XYZ(); @@ -5602,9 +5748,10 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData& data, gp_Ax2d axis( center, vec0 ); gp_Circ2d circ( axis, radius ); - for ( size_t i = iFrom; i < iTo; ++i ) + for ( size_t i = 0; i < _eos.size(); ++i ) { - if ( _eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue; + if ( _eos[i]->Is( _LayerEdge::BLOCKED )) continue; + //if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH )) continue; double newU = uLast * _leParams[i]; gp_Pnt2d newUV = ElCLib::Value( newU, circ ); _eos._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 ); @@ -5639,7 +5786,9 @@ bool _Smoother1D::smoothComplexEdge( _SolidData& data, if ( _offPoints.empty() ) return false; + // ---------------------------------------------- // move _offPoints along normals of _LayerEdge's + // ---------------------------------------------- _LayerEdge* e[2] = { getLEdgeOnV(0), getLEdgeOnV(1) }; if ( e[0]->Is( _LayerEdge::NORMAL_UPDATED )) @@ -5683,7 +5832,9 @@ bool _Smoother1D::smoothComplexEdge( _SolidData& data, } } + // ----------------------------------------------------------------- // project tgt nodes of extreme _LayerEdge's to the offset segments + // ----------------------------------------------------------------- if ( e[0]->Is( _LayerEdge::NORMAL_UPDATED )) _iSeg[0] = 0; if ( e[1]->Is( _LayerEdge::NORMAL_UPDATED )) _iSeg[1] = _offPoints.size()-2; @@ -5747,7 +5898,9 @@ bool _Smoother1D::smoothComplexEdge( _SolidData& data, if ( e[1]->_normal * vDiv1.XYZ() < 0 ) e[1]->_len += d1; else e[1]->_len -= d1; + // --------------------------------------------------------------------------------- // compute normalized length of the offset segments located between the projections + // --------------------------------------------------------------------------------- size_t iSeg = 0, nbSeg = _iSeg[1] - _iSeg[0] + 1; vector< double > len( nbSeg + 1 ); @@ -5771,12 +5924,15 @@ bool _Smoother1D::smoothComplexEdge( _SolidData& data, _offPoints[ _iSeg[0] ]._xyz = pExtreme[0].XYZ(); _offPoints[ _iSeg[1]+ 1]._xyz = pExtreme[1].XYZ(); + // ------------------------------------------------------------- // distribute tgt nodes of _LayerEdge's between the projections + // ------------------------------------------------------------- iSeg = 0; - for ( size_t i = 0; i < _eos._edges.size(); ++i ) + for ( size_t i = 0; i < _eos.size(); ++i ) { - if ( _eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue; + if ( _eos[i]->Is( _LayerEdge::BLOCKED )) continue; + //if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH )) continue; while ( iSeg+2 < len.size() && _leParams[i] > len[ iSeg+1 ] ) iSeg++; double r = ( _leParams[i] - len[ iSeg ]) / ( len[ iSeg+1 ] - len[ iSeg ]); @@ -5785,17 +5941,17 @@ bool _Smoother1D::smoothComplexEdge( _SolidData& data, if ( surface.IsNull() ) { - _eos._edges[i]->_pos.back() = p; + _eos[i]->_pos.back() = p; } else // project a new node position to a FACE { - gp_Pnt2d uv ( _eos._edges[i]->_pos.back().X(), _eos._edges[i]->_pos.back().Y() ); + gp_Pnt2d uv ( _eos[i]->_pos.back().X(), _eos[i]->_pos.back().Y() ); gp_Pnt2d uv2( surface->NextValueOfUV( uv, p, fTol )); p = surface->Value( uv2 ).XYZ(); - _eos._edges[i]->_pos.back().SetCoord( uv2.X(), uv2.Y(), 0 ); + _eos[i]->_pos.back().SetCoord( uv2.X(), uv2.Y(), 0 ); } - SMDS_MeshNode* tgtNode = const_cast( _eos._edges[i]->_nodes.back() ); + SMDS_MeshNode* tgtNode = const_cast( _eos[i]->_nodes.back() ); tgtNode->setXYZ( p.X(), p.Y(), p.Z() ); dumpMove( tgtNode ); } @@ -5836,8 +5992,20 @@ void _Smoother1D::prepare(_SolidData& data) double fullLen = _leParams.back() + pPrev.Distance( SMESH_TNodeXYZ( getLEdgeOnV(1)->_nodes[0])); for ( size_t i = 0; i < _leParams.size()-1; ++i ) _leParams[i] = _leParams[i+1] / fullLen; + _leParams.back() = 1.; } + _LayerEdge* leOnV[2] = { getLEdgeOnV(0), getLEdgeOnV(1) }; + + // get cosin to use in findEdgesToSmooth() + _edgeDir[0] = getEdgeDir( E, leOnV[0]->_nodes[0], data.GetHelper() ); + _edgeDir[1] = getEdgeDir( E, leOnV[1]->_nodes[0], data.GetHelper() ); + _leOnV[0]._cosin = Abs( leOnV[0]->_cosin ); + _leOnV[1]._cosin = Abs( leOnV[1]->_cosin ); + if ( _eos._sWOL.IsNull() ) // 3D + for ( int iEnd = 0; iEnd < 2; ++iEnd ) + _leOnV[iEnd]._cosin = Abs( _edgeDir[iEnd].Normalized() * leOnV[iEnd]->_normal ); + if ( isAnalytic() ) return; @@ -5864,8 +6032,6 @@ void _Smoother1D::prepare(_SolidData& data) _offPoints[i]._param = GCPnts_AbscissaPoint::Length( c3dAdaptor, u0, u ) / _curveLen; } - _LayerEdge* leOnV[2] = { getLEdgeOnV(0), getLEdgeOnV(1) }; - // set _2edges _offPoints [0]._2edges.set( &_leOnV[0], &_leOnV[0], 0.5, 0.5 ); _offPoints.back()._2edges.set( &_leOnV[1], &_leOnV[1], 0.5, 0.5 ); @@ -5903,9 +6069,6 @@ void _Smoother1D::prepare(_SolidData& data) int iLBO = _offPoints.size() - 2; // last but one - _edgeDir[0] = getEdgeDir( E, leOnV[0]->_nodes[0], data.GetHelper() ); - _edgeDir[1] = getEdgeDir( E, leOnV[1]->_nodes[0], data.GetHelper() ); - _leOnV[ 0 ]._normal = getNormalNormal( leOnV[0]->_normal, _edgeDir[0] ); _leOnV[ 1 ]._normal = getNormalNormal( leOnV[1]->_normal, _edgeDir[1] ); _leOnV[ 0 ]._len = 0; @@ -6490,6 +6653,80 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper& } } +//================================================================================ +/*! + * \brief Find _LayerEdge's located on boundary of a convex FACE whose normal + * will be updated at each inflation step + */ +//================================================================================ + +void _ViscousBuilder::findEdgesToUpdateNormalNearConvexFace( _ConvexFace & convFace, + _SolidData& data, + SMESH_MesherHelper& helper ) +{ + const TGeomID convFaceID = getMeshDS()->ShapeToIndex( convFace._face ); + const double preci = BRep_Tool::Tolerance( convFace._face ); + Handle(ShapeAnalysis_Surface) surface = helper.GetSurface( convFace._face ); + + bool edgesToUpdateFound = false; + + map< TGeomID, _EdgesOnShape* >::iterator id2eos = convFace._subIdToEOS.begin(); + for ( ; id2eos != convFace._subIdToEOS.end(); ++id2eos ) + { + _EdgesOnShape& eos = * id2eos->second; + if ( !eos._sWOL.IsNull() ) continue; + if ( !eos._hyp.ToSmooth() ) continue; + for ( size_t i = 0; i < eos._edges.size(); ++i ) + { + _LayerEdge* ledge = eos._edges[ i ]; + if ( ledge->Is( _LayerEdge::UPD_NORMAL_CONV )) continue; // already checked + if ( ledge->Is( _LayerEdge::MULTI_NORMAL )) continue; // not inflatable + + gp_XYZ tgtPos = ( SMESH_NodeXYZ( ledge->_nodes[0] ) + + ledge->_normal * ledge->_lenFactor * ledge->_maxLen ); + + // the normal must be updated if distance from tgtPos to surface is less than + // target thickness + + // find an initial UV for search of a projection of tgtPos to surface + const SMDS_MeshNode* nodeInFace = 0; + SMDS_ElemIteratorPtr fIt = ledge->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face); + while ( fIt->more() && !nodeInFace ) + { + const SMDS_MeshElement* f = fIt->next(); + if ( convFaceID != f->getshapeId() ) continue; + + SMDS_ElemIteratorPtr nIt = f->nodesIterator(); + while ( nIt->more() && !nodeInFace ) + { + const SMDS_MeshElement* n = nIt->next(); + if ( n->getshapeId() == convFaceID ) + nodeInFace = static_cast< const SMDS_MeshNode* >( n ); + } + } + if ( !nodeInFace ) + continue; + gp_XY uv = helper.GetNodeUV( convFace._face, nodeInFace ); + + // projection + surface->NextValueOfUV( uv, tgtPos, preci ); + double dist = surface->Gap(); + if ( dist < 0.95 * ledge->_maxLen ) + { + ledge->Set( _LayerEdge::UPD_NORMAL_CONV ); + if ( !ledge->_curvature ) ledge->_curvature = new _Curvature; + ledge->_curvature->_uv.SetCoord( uv.X(), uv.Y() ); + edgesToUpdateFound = true; + } + } + } + + if ( !convFace._isTooCurved && edgesToUpdateFound ) + { + data._convexFaces.insert( make_pair( convFaceID, convFace )).first->second; + } +} + //================================================================================ /*! * \brief Modify normals of _LayerEdge's on EDGE's to avoid intersection with @@ -6942,6 +7179,8 @@ bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData& data, for ( ; id2face != data._convexFaces.end(); ++id2face ) { _ConvexFace & convFace = (*id2face).second; + convFace._normalsFixedOnBorders = false; // to update at each inflation step + if ( convFace._normalsFixed ) continue; // already fixed if ( convFace.CheckPrisms() ) @@ -7304,6 +7543,59 @@ bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData& data, return true; } +//================================================================================ +/*! + * \brief Return max curvature of a FACE + */ +//================================================================================ + +double _ConvexFace::GetMaxCurvature( _SolidData& data, + _EdgesOnShape& eof, + BRepLProp_SLProps& surfProp, + SMESH_MesherHelper& helper) +{ + double maxCurvature = 0; + + TopoDS_Face F = TopoDS::Face( eof._shape ); + + const int nbTestPnt = 5; + const double oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. ); + SMESH_subMeshIteratorPtr smIt = eof._subMesh->getDependsOnIterator(/*includeSelf=*/true); + while ( smIt->more() ) + { + SMESH_subMesh* sm = smIt->next(); + const TGeomID subID = sm->GetId(); + + // find _LayerEdge's of a sub-shape + _EdgesOnShape* eos; + if (( eos = data.GetShapeEdges( subID ))) + this->_subIdToEOS.insert( make_pair( subID, eos )); + else + continue; + + // check concavity and curvature and limit data._stepSize + const double minCurvature = + 1. / ( eos->_hyp.GetTotalThickness() * ( 1 + theThickToIntersection )); + size_t iStep = Max( 1, eos->_edges.size() / nbTestPnt ); + for ( size_t i = 0; i < eos->_edges.size(); i += iStep ) + { + gp_XY uv = helper.GetNodeUV( F, eos->_edges[ i ]->_nodes[0] ); + surfProp.SetParameters( uv.X(), uv.Y() ); + if ( surfProp.IsCurvatureDefined() ) + { + double curvature = Max( surfProp.MaxCurvature() * oriFactor, + surfProp.MinCurvature() * oriFactor ); + maxCurvature = Max( maxCurvature, curvature ); + + if ( curvature > minCurvature ) + this->_isTooCurved = true; + } + } + } // loop on sub-shapes of the FACE + + return maxCurvature; +} + //================================================================================ /*! * \brief Finds a center of curvature of a surface at a _LayerEdge @@ -7586,13 +7878,14 @@ gp_Ax1 _LayerEdge::LastSegment(double& segLen, _EdgesOnShape& eos) const //================================================================================ /*! - * \brief Return the last position of the target node on a FACE. + * \brief Return the last (or \a which) position of the target node on a FACE. * \param [in] F - the FACE this _LayerEdge is inflated along + * \param [in] which - index of position * \return gp_XY - result UV */ //================================================================================ -gp_XY _LayerEdge::LastUV( const TopoDS_Face& F, _EdgesOnShape& eos ) const +gp_XY _LayerEdge::LastUV( const TopoDS_Face& F, _EdgesOnShape& eos, int which ) const { if ( F.IsSame( eos._sWOL )) // F is my FACE return gp_XY( _pos.back().X(), _pos.back().Y() ); @@ -7601,7 +7894,7 @@ gp_XY _LayerEdge::LastUV( const TopoDS_Face& F, _EdgesOnShape& eos ) const return gp_XY( 1e100, 1e100 ); // _sWOL is EDGE of F; _pos.back().X() is the last U on the EDGE - double f, l, u = _pos.back().X(); + double f, l, u = _pos[ which < 0 ? _pos.size()-1 : which ].X(); Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface( TopoDS::Edge(eos._sWOL), F, f,l); if ( !C2d.IsNull() && f <= u && u <= l ) return C2d->Value( u ).XY(); @@ -7679,7 +7972,7 @@ bool _LayerEdge::SegTriaInter( const gp_Ax1& lastSegment, * \param [in] eov - EOS of the VERTEX * \param [in] eos - EOS of the FACE * \param [in] step - inflation step - * \param [in,out] badSmooEdges - not untangled _LayerEdge's + * \param [in,out] badSmooEdges - tangled _LayerEdge's */ //================================================================================ @@ -9031,6 +9324,11 @@ void _LayerEdge::Block( _SolidData& data ) //if ( Is( BLOCKED )) return; Set( BLOCKED ); + SMESH_Comment msg( "#BLOCK shape="); + msg << data.GetShapeEdges( this )->_shapeID + << ", nodes " << _nodes[0]->GetID() << ", " << _nodes.back()->GetID(); + dumpCmd( msg + " -- BEGIN") + _maxLen = _len; std::queue<_LayerEdge*> queue; queue.push( this ); @@ -9056,6 +9354,8 @@ void _LayerEdge::Block( _SolidData& data ) //if ( edge->_nodes[0]->getshapeId() == neibor->_nodes[0]->getshapeId() ) viscous_layers_00/A3 { newMaxLen *= edge->_lenFactor / neibor->_lenFactor; + // newMaxLen *= Min( edge->_lenFactor / neibor->_lenFactor, + // neibor->_lenFactor / edge->_lenFactor ); } if ( neibor->_maxLen > newMaxLen ) { @@ -9073,6 +9373,7 @@ void _LayerEdge::Block( _SolidData& data ) } } } + dumpCmd( msg + " -- END") } //================================================================================ @@ -9227,6 +9528,7 @@ std::string _LayerEdge::DumpFlags() const case BLOCKED: dump << "BLOCKED"; break; case INTERSECTED: dump << "INTERSECTED"; break; case NORMAL_UPDATED: dump << "NORMAL_UPDATED"; break; + case UPD_NORMAL_CONV: dump << "UPD_NORMAL_CONV"; break; case MARKED: dump << "MARKED"; break; case MULTI_NORMAL: dump << "MULTI_NORMAL"; break; case NEAR_BOUNDARY: dump << "NEAR_BOUNDARY"; break; @@ -9262,7 +9564,7 @@ bool _ViscousBuilder::refine(_SolidData& data) double f,l, u = 0; gp_XY uv; vector< gp_XYZ > pos3D; - bool isOnEdge; + bool isOnEdge, isTooConvexFace = false; TGeomID prevBaseId = -1; TNode2Edge* n2eMap = 0; TNode2Edge::iterator n2e; @@ -9306,6 +9608,9 @@ bool _ViscousBuilder::refine(_SolidData& data) for ( size_t j = 0; j < eos._eosC1[i]->_edges.size(); ++j ) eos._eosC1[i]->_edges[j]->Set( _LayerEdge::SMOOTHED_C1 ); } + isTooConvexFace = false; + if ( _ConvexFace* cf = data.GetConvexFace( eos._shapeID )) + isTooConvexFace = cf->_isTooCurved; } vector< double > segLen; @@ -9321,8 +9626,8 @@ bool _ViscousBuilder::refine(_SolidData& data) if ( eos._sWOL.IsNull() ) { bool useNormal = true; - bool usePos = false; - bool smoothed = false; + bool usePos = false; + bool smoothed = false; double preci = 0.1 * edge._len; if ( eos._toSmooth && edge._pos.size() > 2 ) { @@ -9330,8 +9635,7 @@ bool _ViscousBuilder::refine(_SolidData& data) } if ( smoothed ) { - if ( !surface.IsNull() && - !data._convexFaces.count( eos._shapeID )) // edge smoothed on FACE + if ( !surface.IsNull() && !isTooConvexFace ) // edge smoothed on FACE { useNormal = usePos = false; gp_Pnt2d uv = helper.GetNodeUV( geomFace, edge._nodes[0] ); -- 2.30.2