[wyldckat]

Hi Sandeep,

Many thanks for the feedback! I should have posted the several tests I've made so far, because I did try several different ways of mapping the "pointDisplacement" field and respective mesh.

That one in particular I think didn't work, although it's possible that I didn't manage to define the proper usage of each object and respective methods.

Later today I'll post the attempts I've made at mapping the field.

In the mean time, the problem is also hat I did further (re)searching on this topic, which lead me to thinking that I would need several more complex hacks that I found in OpenFOAM's source code, specifically the ones in:

• "fvMeshDistribute.C":

  00544     // Store boundary fields (we only do this for surfaceFields)
  00545     PtrList<FieldField<fvsPatchField, scalar> > sFlds;
  00546     saveBoundaryFields<scalar, surfaceMesh>(sFlds);
  ...
  00555 
  00556     // Change the mesh (no inflation). Note: parallel comms allowed.
  00557     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
  00558 
  00559     // Update fields. No inflation, parallel sync.
  00560     mesh_.updateMesh(map);
  00561 
  00562     // Map patch fields using stored boundary fields. Note: assumes order
  00563     // of fields has not changed in object registry!
  00564     mapBoundaryFields<scalar, surfaceMesh>(map, sFlds);
  ...

  More specifically, I would need to borrow code from "fvMeshDistributeTemplates.C" and "fvMeshDistribute.C", for specifically mapping the point mesh objects, instead of volume and surface meshes.
• Then there is:
  

  void Foam::displacementFvMotionSolver::updateMesh(const mapPolyMesh& mpm)

  But seems even more complex and only details about point movement...
• And I know there is a similar hack to the one you've showed, but made for surface meshes... but I can't find it.

Like I've said, later today I'll post the attempts I've made so far.


Best regards,
Bruno

[wyldckat]

Hi Sandeep,

OK, I've got mixed results...

First for the partially-good news:

While using the hack to use the file system and by adding the sphere (named "oggetto") to the "noModificationPatches" and "noSwapPatches", it does allow to keep the sphere's surface in perfect condition while it moves.
The downside is that it crashes when the cells in the way start to get too tight and eventually collapse into negative volumes.

The error message:

Courant Number mean: 0.04541765902 max: 4.191785362
deltaT = 3.333333333e-05
Time = 0.0202

sixDoFRigidBodyMotion constraints converged in 1 iterations
Constraint force: (0.6537136637 0.02252467437 0)
Constraint moment: (0 0 0)
Centre of mass: (-1.842749434e-08 1.08451596e-08 0.2001177404)
Linear velocity: (-0.0004352669538 8.108849543e-05 0.04280597798)
Angular velocity: (1.573833728 -1.569414796 -0.03449086776)

************** QualityAssessor(free only) Summary **************

  Evaluating quality for 72216 free elements of 72375 total elements.
  THERE ARE 1 INVERTED ELEMENTS. 
  (Elements invalid at 1 of 72216 sample locations.)

  1 OF 72216 ENTITIES EVALUATED TO AN UNDEFINED VALUE FOR AspectRatioGamma

          metric     minimum     average         rms     maximum    std.dev.
AspectRatioGamma 1.002851198   15.115086 3721.202667     1000000 3721.171969

 Warning: Minimum cell quality is: -0.07980588317 at cell: 68575

 ~~~ Mesh Quality Statistics ~~~ 
 Min: -0.07980588317
 Max: 0.9981037065
 Mean: 0.8637096283
 Cells: 72375
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 



--> FOAM FATAL ERROR: 
Encountered negative cell-quality!
Edge: 89295: (3253 3261)
vertexHull: 5(227 9861 10215 9856 3297)
Minimum Quality: -0.07980588317

    From function scalar dynamicTopoFvMesh::computeMinQuality(const label eIndex, labelList& hullVertices) const
    in file edgeSwap.C at line 2246.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) in "/home/bmss/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#1  Foam::error::abort() in "/home/bmss/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2  Foam::dynamicTopoFvMesh::computeMinQuality(int, Foam::List<int>&) const in "/home/bmss/OpenFOAM/bmss-2.1.x/platforms/linux64GccDPOpt/lib/libdynamicTopoFvMesh.so"
#3  Foam::dynamicTopoFvMesh::swap3DEdges(void*) in "/home/bmss/OpenFOAM/bmss-2.1.x/platforms/linux64GccDPOpt/lib/libdynamicTopoFvMesh.so"
#4  Foam::dynamicTopoFvMesh::threadedTopoModifier() in "/home/bmss/OpenFOAM/bmss-2.1.x/platforms/linux64GccDPOpt/lib/libdynamicTopoFvMesh.so"
#5  Foam::dynamicTopoFvMesh::update() in "/home/bmss/OpenFOAM/bmss-2.1.x/platforms/linux64GccDPOpt/lib/libdynamicTopoFvMesh.so"
#6  
 in "/home/bmss/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/bin/pimpleDyMFoam"
#7  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#8  
 in "/home/bmss/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/bin/pimpleDyMFoam"
Aborted (core dumped)

The "dynamicMeshDict":

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM Extend Project: Open Source CFD        |
|  \\    /   O peration     | Version:  1.6-ext                               |
|   \\  /    A nd           | Web:      www.extend-project.de                 |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      dynamicMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

//- Select the type of dynamicFvMesh
dynamicFvMesh           dynamicTopoFvMesh;

//- Select the type of motionSolver
motionSolverLibs ("libmesquiteMotionSolver.so" "libfvMotionSolvers.so");
solver                  mesquiteMotionSolver;

mesquiteOptions
{
    usePointDisplacement yes;
  
    //- Optimization metric
    optMetric               AspectRatioGamma;

    //- Objective function
    objFunction             LPtoP;

    //- Optimization algorithm
    optAlgorithm            FeasibleNewton;

    //- Termination criteria sub-dictionaries
    //- (takes default values if not specified)
    //-   Specifying an empty sub-dictionary
    //-   terminates with available options
    tcInner
    {
        absGradL2            1e-4;
        cpuTime              0.5;
        // lele: iterationLimit 10;
    }

    // tcOuter
    // {}

    //- For composite functions, two objectives need to be specified
    // firstFunction           LPtoP;
    // secondFunction          LInf;

    //- For scaled functions, scale and objective needs to be specified
    // scaleFunction        PMeanP;
    // scale                1.5;

    //- Power value for the LPtoP objective function
    pValue                  2;
    power                   2;

    //- Specify a tolerance for the surface-smoothing CG solver
    tolerance               1e-2;

    //- Specify number of CG sweeps for surface-smoothing
    nSweeps                 2;

    //- Specify slip patches for the motionSolver
    slipPatches
    {
        //wall;
        //oggetto;
    }

    //- Constrain surface mesh-motion on a specified cylinder
    cylindricalConstraints
    {
        //- Specify options per slip patch
    //    sideWall
    //    {
    //        axisPoint      (0.0 0.0 0.0);
    //        axisVector     (0.0 0.0 1.0);
    //        radius          1.0;
    //    }
    }

    //- Specify fixedValue patches for the motionSolver
    fixedValuePatches
    {
      oggetto
      {
          type            sixDoFRigidBodyDisplacement;
          centreOfMass    (0 0 0.2);
          momentOfInertia (0.000001 0.000001 0.000001);
          mass            0.01;
          rhoName         rhoInf;
          rhoInf          1000;  // needed only for solvers solving for kinematic pressure
          report          on;
          value           uniform (0 0 0);
          constraints
          {
            maxIterations 500000000;
            fixedLine1
            {
              sixDoFRigidBodyMotionConstraint fixedLine;
              tolerance 1e-6;
              relaxationFactor 0.7;
              fixedLineCoeffs
              {
                refPoint (0 0 0.2);
                direction (0 0 1);
              }
            }
          }
      }
    }
    
    //Other fixedValuePatch types can be found here:
    //src/dynamicMesh/meshMotion/fvMotionSolver/pointPatchFields/derived
    
    // angularOscillatingDisplacement
    // angularOscillatingVelocity
    // oscillatingDisplacement
    // oscillatingVelocity
    // surfaceDisplacement
    // surfaceSlipDisplacement

    //- Specify interval for surface smoothing
    surfInterval            1;
}

//- Options for dynamicTopoFvMesh
dynamicTopoFvMesh
{
    //- Should all options be made mandatory?
    //- Useful for first-time use.
    allOptionsMandatory no;

    //- Set run-time debug level [0-5]
    //debug               2;

    //- Specify the number of threads
    threads             1;

    //- Specify re-meshing interval
    //- Negative value implies no re-meshing
    interval            1;

    //- Specify whether the length-scale field
    //- should be dumped to disk
    dumpLengthScale     false;

    //- sliverThreshold specifies the
    //- quality criteria for sliver removal.
    sliverThreshold     0.35;

    //- Should the tool attempt to remove slivers
    //- that fall below the sliverThreshold value?
    removeSlivers       true;

    //- Skip mapping step. Useful while using
    //- this tool as a pre-processor
    // skipMapping         true;

    // Toggle edgeRefinement on/off
    edgeRefinement      no;

    //- Options for edge-bisection/collapse.
    //-   The edgeRefinement flag must be set for
    //-   the following options to have effect
    refinementOptions
    {
        collapseRatio   0.5;
        bisectionRatio  1.5;
        growthFactor    1.0;

        //- By default, existing boundary edge-lengths
        //- are used for length-scales.
        //- Length-scale can be fixed for certain patches.
        fixedLengthScalePatches
        {
            oggetto     0.2;
        }

        //- Avoid refinement on certain patches, if necessary
        noModificationPatches
        {
            inlet;
            outlet;
            wall;
            oggetto;
        }

        FreeSurfacePatch
        {
            oggetto;
        }
        
        
        //- Set floating length-scale values on certain patches
        freeLengthScalePatches
        {
            inlet;
            outlet;
            wall;
        }

        //- Limit lengthScales to specified values, if necessary
        // minLengthScale   0.1;
        // maxLengthScale   0.3;

        //- Field-based refinement options
        // fieldRefinement  gamma;
        // fieldLengthScale 0.005;
        // lowerRefineLevel 0.001;
        // upperRefineLevel 0.999;
        // maxRefineLevel   4;
        // meanScale        0.015;
    }

    //- If the number of modifications are to be limited, set this option
    // maxModifications   1000;

    //- Load custom libraries for metrics
    // tetMetricLibs      ("libtetMetrics.so");

    //- Tetrahedral mesh quality metric
    tetMetric          Knupp;

    //- Avoid 2-2 swapping on certain patches
    noSwapPatches
    {
      oggetto;
    }
}

// ************************************************************************* //

Attached are the two images that show the initial position and the final position nearly before the crash, due to a collapsed cell.

The bad news:

I still haven't been able to make the proper update to work. I've tried to use the following code, but without success (this is the patch for the code):

----------------- mesquiteMotionSolver/mesquiteMotionSolver.C -----------------
index 468be60..a482416 100644
@@ -4738,6 +4738,18 @@ void mesquiteMotionSolver::update(const mapPolyMesh& mpm)
         Info<< "Clearing out mesquiteMotionSolver for topo-changes" << endl;
     }
 
+    // pointMesh
+    if (Mesh_.thisDb().foundObject<pointMesh>(pointMesh::typeName))
+    {
+        const_cast<pointMesh&>
+        (
+            Mesh_.thisDb().lookupObject<pointMesh>
+            (
+                pointMesh::typeName
+            )
+        ).updateMesh(mpm);
+    }
+
     motionSolver::updateMesh(mpm);
 
     if (surfaceSmoothing_)

The crash occurs just the same, namely whenever the "boundaryConditions_()" are accessed after the mesh has been changed.

Other attempts I did in the past were, both in "mesquiteMotionSolver::update(...)":

• Force points to move:

      const_cast<pointMesh&>(boundaryConditions_().mesh())
          .movePoints(Mesh_.points());

• Force an update of the mesh:

      const_cast<pointMesh&>(boundaryConditions_().mesh())
          .updateMesh(mpm);

Both crashed while trying to make the move, presumably due to a non-existing mesh to be moved...

Best regards,
Bruno

[deepsterblue]

Bruno,

You could have the smoother work a little harder at each time step, and see if the simulation goes further. Specify a higher cpuTime value for tcInner and see if that helps. Also, you have edgeRefinement set to 'no', which means that cells cannot be collapsed, and that probably explains why you see negative volumes.

You can play around with the growth factor to coarsen the mesh away from the 'oggetto' patch (not too agressive, otherwise you'll get too many collapses, which you can restrict using the maxModifications entry).

I'll look into the mapping when I have some time on my hands.

Sandeep

[wyldckat]

Hi Sandeep,

Many thanks for the hints!

Unfortunately, by removing the cells that were in the way, the saved "pointDisplacement" of the internal field has a different number of cells, which leads the solver stopping to complain about the different array sizes.

Nonetheless, this does give me (again) the idea for adopting a hack similar to the strategy used for "fixedValuePatches"...

Best regards,
Bruno

[deepsterblue]

Bruno,

Can you give me a test case that I can use? Also, you'll have to list out the details of how to reproduce the problem.

Sandeep

[wyldckat]

Hi Sandeep,

Many thanks! Although I'll only be able to provide a clean/smaller case and more explicit instructions in about 5 hours from now


In the mean time, the case is the one mentioned in post #13, namely: https://www.dropbox.com/s/bs3ehwn6xq...topo_v2.tar.gz

Unpack and simply run pimpleDyMFoam.
It should crash as soon as "boundaryConditions_().correctBoundaryConditions() " is called the second time around, if I'm not mistaken.
moveDynamicMesh is not enough, due to the need for forces to work, which depend on "U" and "p"...


This is assuming that the whole "dynamicTopoFvMesh" code only has the first two source code patches you provided on this thread, for using with OpenFOAM 2.1.x.

Unless you prefer the hack I tried to use, relying on the file system, for which the "pointDisplacement" file was placed in the "constant" folder; the source code patch is on post #19.

Best regards,
Bruno

[deepsterblue]

Bruno,

I've committed a fix and your case seems to be running with point fields mapped. I haven't tried extensively with edgeRefinement on, so you might want to test it out.

Unfortunately, since your length scales vary quite drastically across the 'wall' patch, the algorithm has a hard time trying to figure out what the optimal length scale is, and so you tend to get excessive bisections / collapses. You have a length scale for 'oggetto' specified as 0.2, while the paraview ruler tells me it's 0.002 - so you may want to check that.

Until I have a more general solution, one workaround is to limit the number at each time-step using the maxModifications entry, or perhaps you could split up the patch into several smaller ones with a uniform length scale for each.

Thanks for the bug-report.

[wyldckat]

Hi Sandeep,

You're very welcome about the bug report! And many thanks for the fix and for the suggestions!

And no wonder I couldn't figure it out... I was only looking at the source code for "mesquiteMotionSolver"

Best regards,
Bruno

[joegi.geo]

Hi Sandeep,

I wonder if there is a difference between the implementation of patch motion between dynamicTopoFvMesh and dynamicMotionSolverFvMesh?. I am trying to do a simple oscillating body simulation by using both libraries and the output of the patch motion is totally different (about an order of magnitude in amplitude).

Using dynamicMotionSolverFvMesh i set the moving patch in pointDisplacement as follows:

    wing
    {
        type            oscillatingDisplacement;
        value           uniform ( 0 0 0 );
        amplitude       ( 0 0.2 0 );
        omega           6.28318;
    }

which gives me the desired output.

Using dynamicTopoFvMesh/mesquite I set the moving patch in dynamicMeshDict exactly in the same way, by just adding

    wing
    {
        type            oscillatingDisplacement;
        value           uniform ( 0 0 0 );
        amplitude       ( 0 0.2 0 );
        omega           6.28318;
    }

to dynamicMeshDict in the mesquiteOptions entry,and the results are totally different. I also tried by using the pointDisplacement dictionary and the output didn't change.

Any clue?,

jg

[deepsterblue]

I see that both your entries are the same, but how different is "totally different"? What was your output? If you post your entire dynamicMeshDict file, that would be helpful.

[joegi.geo]

Hi,

Ok, this is the entry where i set the patch motion

solver                  mesquiteMotionSolver;

dynamicFvMesh dynamicTopoFvMesh;

dynamicFvMeshLibs         ("libdynamicTopoFvMesh.so");

motionSolverLibs  ("libfvMotionSolvers.so"  "libmesquiteMotionSolver.so");

mesquiteOptions
{
    // Optimization metric
    optMetric               AspectRatioGamma;

    // Objective function
    objFunction             LPtoP;

    // Optimization algorithm
    optAlgorithm            FeasibleNewton;

    // Termination criteria sub-dictionary (takes default values if not specified)
    // Specifying an empty sub-dictionary terminates with available options
    tcInner
    {
       absGradL2            1e-4;
       cpuTime              1.25;
    }
    // tcOuter
    // {}

    // For composite functions, two objectives need to be specified
    firstFunction           LPtoP;
    secondFunction          LInf;

    // For scaled functions, scale and objective needs to be specified
    // scaleFunction        PMeanP;
    // scale                1.5;

    // Power value for the LPtoP objective function
    pValue                  2;
    power                   2;

    // Specify a tolerance for the CG solver
    tolerance               1e-2;

    // Specify number of CG sweeps
    nSweeps                 1;

    // Specify slip patches for the motionSolver
    slipPatches
    {

    }

    // Set run-time debug level
    debug                   0;

    // Specify interval for surface smoothing
    surfInterval            1;



    //- Specify fixedValue patches for the motionSolver
    fixedValuePatches
    {

        ball
        {
        type            oscillatingDisplacement;
            amplitude     (0.1 0 0);
            omega         6.28316;
            value         uniform (0 0 0);
    }

    }

}

If dynamicTopoFvMesh uses the same patch motion class as dynamicMotionSolverFvMesh (which I think it does), I should get the same result for both dynamic meshing libraries, am i right?. In my case the output of both classes are totally different.

Even if I set the entry usePointDisplacement yes, and use the pointDisplacement dictionary (same one as the one used with dynamicMotionSolverFvMesh), the patch motion of both libraries are different.

[deepsterblue]

And what's the output?

[joegi.geo]

The simulation runs fine (I am using moveDynamicMesh), the smoothing works fine, the remeshing works fine. But, nevertheless the fact that I am using the same patch motion for both dynamic meshing classes (libdynamicTopoFvMesh and dynamicMotionSolverFvMesh), I think the motion I am getting with libdynamicTopoFvMesh is not the right one, the amplitude is like one order of magnitude higher compare to the output of dynamicMotionSolverFvMesh.

[deepsterblue]

Ah... I know what the problem is:

It appears that the oscillatingDisplacement patchField provides an absolute value of displacement (from the mean), while the applyFixedValuePatches() function in mesquiteMotionSolver expects an incremental displacement vector for each time-step.

You see the acceleration because the displacement is continuously updated as:

refPoints_ += dPointField; in mesquiteMotionSolver.C

You could write a new boundary condition similar to oscillatingDisplacement which provides displacement increments, or alternatively, I also provide a timeVaryingDisplacement patchField in the repository that uses an interpolation table from file, which might be easier. I would have to look into a fix for the absolute vs. incremental confusion.

Can you try changing the "+=" operator in applyFixedValuePatches to simply "=" and see if that makes a difference?

[joegi.geo]

It crash immediately.

Let's see if I can write the bc, in anycase I am not in a hurry, so probably you will come with a solution first.

Does this affect the behavior of the 6DOF solver? it is safe to use?


Have a nice day,

joel

[deepsterblue]

Can you post a test case with dynamicMeshDict entries for both solver types?

[joegi.geo]

Hi Sandeep,

You can use the ballTranslation tutorial from the last openfoam workshop. Use the attached files to run the case using dynamicMotionSolverFvMesh.
I tested everything and runs fine.

Running this two cases you clearly see the difference in the patch motion between the two solvers.

Have a nice day,

joel

[deepsterblue]

joel,

Can you try the angularOscillating variant with zero angular velocity and see if there's a difference between both solvers?

Thanks

[joegi.geo]

Hi Sandeep,

Same problem with angularOscillatingDisplacement.


Joel

[joegi.geo]

Hi Sandeep,

Any update regarding this issue?