Is there a way to do volume mesh deformation using runScript?

[Mohsin-Bashir]

In the runScript, I'm looking for a way to execute volume mesh deformation only if I can update the design surface mesh coordinates using my own module.

I'm aware that pyGeo supplies design surface mesh coordinates to Idwarp for volume mesh deformation. However, I aim to bypass the use of pyGeo and directly manage Idwarp through the runScript. e.g. for running primal only, the following function is executed:
optFuncs.runAdjoint()
Is there such a function for volume mesh deformation too?
Any help regarding this will be much appreciated.

[friedenhe]

This line deforms the FFD point and returns the deformed surface mesh coordinate

dafoam/dafoam/pyDAFoam.py

Line 940 in 0f42fe5

xs = self.DVGeo.update(self.ptSetName, config=None)

Then, it gives the deformed surface point to DAFoam by calling line 947. Finally, it calls line 952 to deform the volume mesh.

You can replace line 940 with your own geometry module.

[Mohsin-Bashir]

This requires changing the source code. However, I am looking for a way to do it without changing source code.

I wonder if I can use a script similar to runScript_v2.py to do volume mesh deformation.
Here is my script deform.py deform.py.txt
I used following extra lines (in script) to do volume mesh deformation but it did not change the points.gz file.

meshSurfCoords=DASolver.getSurfaceCoordinates()
DASolver.setSurfaceCoordinates(meshSurfCoords)
DASolver.mesh.warpMesh()

Workflow:
My geometry module outputs coordinates of deformed design surface in a file. I used that file to update points.gz file (changed unreformed design surface coordinates with deformed ones). After this, I ran deform.py but it did not change points.gz file, though ran without any issues deform_py output.txt. (I was expecting that running deform.py will change the coordinates in points.gz file by doing volume deformation)

Please guide me if its possible to do this without changing source code.

[Mohsin-Bashir]

I guess there is a function named writeGrid for this in idwarp code as shown below: https://github.com/mdolab/idwarp/blob/main/idwarp/UnstructuredMesh.py#L479

But when I use this in my script
DASolver.mesh.writeGrid()
It gives me this error:

Traceback (most recent call last):
  File "deform.py", line 203, in <module>
    DASolver.mesh.writeGrid()
  File "/home/dafoamuser/dafoam/packages/miniconda3/lib/python3.8/site-packages/idwarp/UnstructuredMesh.py", line 505, in writeGrid
    self._writeOpenFOAMVolumePoints(self.getCommonGrid())
AttributeError: 'USMesh' object has no attribute '_writeOpenFOAMVolumePoints'

[friedenhe]

No, modifying the points.gz file won't work because DAFoam does not read points.gz during optimization. It only writes points.gz. In other words, the mesh information is stored in memory during optimization, so any changes to the points.gz will not impact how DAFoam deforms the mesh. In this sense, you have to change the source code. You need to convert your surface mesh to a format that is consistent with meshSurfCoords, and then provide it to the Python layer (the line you cited above).

[Mohsin-Bashir]

I see.
And after changing source code, how can I just run volume mesh deformation? Because at a time I just want to run volume deformation not the whole optimization.

[friedenhe]

You can run the flow for only one step. It will output the deformed points.gz in the output folder/polyMesh.

[Mohsin-Bashir]

This line deforms the FFD point and returns the deformed surface mesh coordinate

dafoam/dafoam/pyDAFoam.py

Line 940 in 0f42fe5

xs = self.DVGeo.update(self.ptSetName, config=None)

Then, it gives the deformed surface point to DAFoam by calling line 947. Finally, it calls line 952 to deform the volume mesh.

You can replace line 940 with your own geometry module.

Hi I have made some changes to the source code (specifically pyDAFoam.py file).
xs = self.DVGeo.update(self.ptSetName, config=None)
I changed the above line to xs=self.getwarpedpoints() and added this new function to pyDAFoam which is as follows:

    def getwarpedpoints(self):
        #reads data from surface_points.data file. Divides each coordinate by 1000 and then convert it into numpy array
        # File path
        surf_filename = '/home/dafoamuser/mount/oneraM6/surface_positions.dat'
        
                
        data = []
        with open(surf_filename, 'r') as file:
            for line in file:
                parts = line.strip().split()
                if len(parts) == 4:
                    index = int(parts[0])
                    x = float(parts[1]) / 1000
                    y = float(parts[2]) / 1000
                    z = float(parts[3]) / 1000
        
                    data.append([index, x, y, z])        
        return np.array(data)[:,1:]

Now the problem is, when I run my case in parrallel, xs=self.getwarpedpoints() gives the whole array all at once instead of first dividing the array among different processors as it was doing earlier for xs = self.DVGeo.update(self.ptSetName, config=None)
Could you please guide me how to deal with this problem?
pyDAFoam.py.txt
logOpt.txt

[friedenhe]

If you use your own surface mesh deformation function, you have to manually split the surface points into multiple processors, following OpenFOAM's standard. xs should be the surface point coordinates owned by local processors

[Mohsin-Bashir]

I don't know what pattern OpenFOAM follows to distribute surface point coordinates among multiple processors. Could you please link some relevant documentation or other resources to understand this?
Currently, I have my surface point coordinates in this format: surface_positions.txt
Update: I have found a way to split surface points among multiple processors:

1. I used the nodeIDs from this file processorN/constant/polyMesh/pointsProcAddressing.gz. But this contains all the node IDs not just related to design surface.
2. So, I just took the intersection of all node ids on a processor with the all node ids. (to get node ids related to only design surface on a processor)
3. Finally, I read the corresponding coordinates from my surface_positions.dat file for each processor.

I don't know how efficient this method is, so, still looking for other methods.

[friedenhe]

This is more like an OpenFOAM question. The way you map the local (owned by each processor in parallel) and global (serial runs) node IDs is correct. To figure out which node ID belongs to the design surface, you have to use label pointIndex = mesh.boundaryMesh()[[patchIndex][faceIndex][pointIndex]; This will give you the node ID for a specific surface node determined by the patchIndex, faceIndex, and pointIndex.

[Mohsin-Bashir]

Thank you Professor!!