Virtual Element Method #1220

andrea-borio · 2020-10-27T16:23:03Z

First hypothetical structure for VEM classes.

Also added a test to check build.

…ementBase.cpp The generation of the test cube is now done succefully.

…tors

…rojectors. Add computation of rotated points and diameter. Add calls to get the required geometrical quantities.

…rojectors. Add computation of integrals of 2D monomials on the boundary of the face.

The method now fills a vector containing the integral mean of projections of basis functions whose support intersects the face.

Bug-fixes and changes to face computations.

Add computation of integral mean of derivatives of basis functions. Add computation of the stabilization matrix. Add computation of the number of quadrature points. Code cleanups.

Fixed bugfix in computation of basisTimesMonomNormalDerBoundaryInt in ConformingVirtualElementOrder1::ComputeProjectors(). Changed inputs of ComputeFaceIntegrals() consequently. All properties of ConformingVirtualElementOrder1 set to public to allow easier testing. To be changed later, when getters are defined.

Face rotation matrices are now computed on-the-fly using pre-computed face normals.

andrea-borio · 2020-12-10T15:54:02Z

This is the document with the pseudo-code of the implementation.

main.pdf

Using arraySlices caused build errors with CUDA.

francoishamon

Thank you very much @andrea-borio , the PR looks great! I just made some small comments, but most of them can be addressed in the next PR if that's easier.

francoishamon · 2021-01-11T19:13:30Z

src/coreComponents/virtualElement/ConformingVirtualElementOrder1.hpp

+  localIndex numQuadraturePoints;
+  localIndex numSupportPoints;
+  array1d< real64 > basisFunctionsIntegralMean;
+  array2d< real64 > basisDerivativesIntegralMean;
+  array2d< real64 > stabilizationMatrix;


Suggested change

localIndex numQuadraturePoints;

localIndex numSupportPoints;

array1d< real64 > basisFunctionsIntegralMean;

array2d< real64 > basisDerivativesIntegralMean;

array2d< real64 > stabilizationMatrix;

localIndex m_numQuadraturePoints;

localIndex m_numSupportPoints;

array1d< real64 > m_basisFunctionsIntegralMean;

array2d< real64 > m_basisDerivativesIntegralMean;

array2d< real64 > m_stabilizationMatrix;

Do they need to be public? Can they be protected instead?

Thank you François! These properties should be protected for sure. They are public now for testing and because we still need to agree on suitable getters.

src/coreComponents/virtualElement/unitTests/CMakeLists.txt

src/coreComponents/virtualElement/VirtualElementBase.hpp

francoishamon · 2021-01-11T19:59:19Z

src/coreComponents/virtualElement/ConformingVirtualElementOrder1.cpp

+
+  // Compute other geometrical properties.
+  //  - compute rotation matrix.
+  array2d< real64 > faceRotationMatrix( 3, 3 );


Suggested change

array2d< real64 > faceRotationMatrix( 3, 3 );

real64 faceRotationMatrix[ 3 ][ 3 ] = {{ 0.0 }};

francoishamon · 2021-01-11T20:10:22Z

src/coreComponents/virtualElement/ConformingVirtualElementOrder1.cpp

+  array1d< real64 > monomInternalIntegrals( 2 );
+  threeDMonomialIntegrals.resize( 3 );
+  LvArray::tensorOps::fill< 2 >( monomInternalIntegrals, 0.0 );
+  LvArray::tensorOps::fill< 3 >( threeDMonomialIntegrals, 0.0 );


Here to avoid resizing and filling the vectors with zero you could just do something more concise I think:

real64 monomInternalIntegrals[ 2 ] = { 0.0 };

and define in the other function threeDMonomialIntegrals as

real64 threeDMonomialIntegrals[ 3 ] = { 0.0 };

francoishamon · 2021-01-11T21:02:41Z

src/coreComponents/virtualElement/ConformingVirtualElementOrder1.cpp

+  // - compute cell diameter.
+  real64 cellDiameter = 0;
+  for( localIndex numVertex = 0; numVertex < numCellPoints; ++numVertex )
+  {
+    for( localIndex numOthVertex = 0; numOthVertex < numVertex; ++numOthVertex )
+    {
+      array1d< real64 > vertDiff( 3 );
+      LvArray::tensorOps::copy< 3 >( vertDiff, nodesCoords[cellToNodes( cellIndex, numVertex )] );
+      LvArray::tensorOps::subtract< 3 >( vertDiff,
+                                         nodesCoords[cellToNodes( cellIndex, numOthVertex )] );
+      real64 const candidateDiameter = LvArray::tensorOps::l2NormSquared< 3 >( vertDiff );
+      if( cellDiameter < candidateDiameter )
+        cellDiameter = candidateDiameter;
+    }
+  }
+  cellDiameter = LvArray::math::sqrt< real64 >( cellDiameter );
+  real64 const invCellDiameter = 1.0/cellDiameter;


I suggest that in the next PR we move this piece of code computing the cell diameter could be moved to ComputationalGeometry since it seems generic and could be useful for others.

francoishamon · 2021-01-11T21:06:09Z

src/coreComponents/virtualElement/ConformingVirtualElementOrder1.cpp

+  array1d< real64 > monomInternalIntegrals( 3 );
+  LvArray::tensorOps::fill< 3 >( monomInternalIntegrals, 0.0 );


These two lines could be simplified to:

real64 monomInternalIntegrals[ 3 ] = { 0.0 };

src/coreComponents/virtualElement/ConformingVirtualElementOrder1.cpp

francoishamon · 2021-01-11T21:46:58Z

src/coreComponents/virtualElement/ConformingVirtualElementOrder1.cpp

+void ConformingVirtualElementOrder1::ComputeProjectors( MeshLevel const & mesh,
+                                                        localIndex const & regionIndex,
+                                                        localIndex const & subRegionIndex,
+                                                        localIndex const & cellIndex )


In the next PR when you want to call this function inside a kernel to compute the basis functions and the stabilization matrix in each cell, you may have to update the signature of this function a little bit.

Instead of passing the mesh and the triplet of indices that identifies a cell, I think you would directly pass the volume and center of the particular cell, plus some arrayViews (say, faceCenters, faceNormals, edgeToNodes, etc), arraySlices (say, elemToFaces and cellToNodes for the particular cell) that you need inside the function. You would also pass a reference to the three arrays (basis functions, derivatives, and stabilization matrix) that you want to populate in this function. With these changes (and maybe some small others) it will be possible to call ComputeProjectors inside a kernel.

Yes, I agree. We will make use of a StackVariables struct, which is set up in a kernel launch looping over cells of a subregion before entering the quadrature point kernel (see kernelLaunch in the base interface for finite element kernels). For example, in SolidMechanicsSmallStrainQuasiStaticKernel.hpp:

stack arrays for the vertex coordinates, displacement, ... are added to StackVariables (see)

values are copied in the setup method (see)

We need to carefully evaluate what we should copy and what we should recompute.

francoishamon · 2021-01-11T21:50:22Z

src/coreComponents/virtualElement/ConformingVirtualElementOrder1.cpp

+  array1d< real64 > basisBoundaryIntegrals( numCellPoints );
+  array2d< real64 > basisTimesNormalBoundaryInt( 3, numCellPoints );
+  array2d< real64 > basisTimesMonomNormalDerBoundaryInt( numCellPoints, 3 );
+  array1d< real64 > monomBoundaryIntegrals( 4 );


Suggested change

array1d< real64 > monomBoundaryIntegrals( 4 );

real64 monomBoundaryIntegrals[ 4 ] = { 0.0 };

and in general I think you can replace (maybe in the next PR) all the array1d and array2d that have a fixed size with C-arrays.

castelletto1

Thanks @andrea-borio! The PDF with the pseudo-code of the implementation was extremely useful. Let's try to have this PR merged as soon as possible and start the next one right after.

If the usual unit test keeps failing for Pangea2, we might have to deactivate them again (after consulting with @TotoGaz).

src/coreComponents/virtualElement/ConformingVirtualElementOrder1.cpp

castelletto1 · 2021-01-15T17:46:58Z

src/coreComponents/virtualElement/ConformingVirtualElementOrder1.hpp

+                             array1d< real64 > & threeDMonomialIntegrals );
+
+public:
+  ConformingVirtualElementOrder1() {}


Suggested change

ConformingVirtualElementOrder1() {}

My compiler says constructors cannot be declared ‘virtual’ [-fpermissive]

Correct, remove the constructor all together.

src/coreComponents/virtualElement/ConformingVirtualElementOrder1.hpp

castelletto1 · 2021-01-15T19:24:45Z

src/coreComponents/virtualElement/ConformingVirtualElementOrder1.cpp

+void ConformingVirtualElementOrder1::ComputeProjectors( MeshLevel const & mesh,
+                                                        localIndex const & regionIndex,
+                                                        localIndex const & subRegionIndex,
+                                                        localIndex const & cellIndex )


Yes, I agree. We will make use of a StackVariables struct, which is set up in a kernel launch looping over cells of a subregion before entering the quadrature point kernel (see kernelLaunch in the base interface for finite element kernels). For example, in SolidMechanicsSmallStrainQuasiStaticKernel.hpp:

stack arrays for the vertex coordinates, displacement, ... are added to StackVariables (see)

values are copied in the setup method (see)

We need to carefully evaluate what we should copy and what we should recompute.

TotoGaz · 2021-01-16T02:09:09Z

If the usual unit test keeps failing for Pangea2, we might have to deactivate them again (after consulting with @TotoGaz).

Yes, deactivate what you need. Just open an issue so we keep track of it.

…tify_style.

andrea-borio · 2021-01-19T16:16:11Z

I labeled the PR as ready to be merged even though there are some issues with Travis. Let's discuss this at today's meeting.

andrea-borio added 21 commits October 27, 2020 17:03

First hypothetical structure for VEM classes.

6dbee70

Also added a test to check build.

Merge remote-tracking branch 'origin/develop' into feature/borio/vem

55a55d5

Modify parameters of VirtualElementBase::ComputeProjectors()

40548a0

Trying to get a test mesh in testVirtualElementBase.cpp

29d38c8

Add main to src/coreComponents/virtualElement/unitTests/testVirtualEl…

aaaac41

…ementBase.cpp The generation of the test cube is now done succefully.

Merge remote-tracking branch 'origin/develop' into feature/borio/vem

7eecf73

Start implementation of ConformingVirtualElementOrder1::ComputeProjec…

0a335ab

…tors

Merge remote-tracking branch 'origin/develop' into feature/borio/vem

2ee134f

WIP on implementation of ConformingVirtualElementOrder1::ComputeFaceP…

cb770c3

…rojectors. Add computation of rotated points and diameter. Add calls to get the required geometrical quantities.

Bugfix in computation of face diameter.

0e739b6

WIP on implementation of ConformingVirtualElementOrder1::ComputeFaceP…

8f13824

…rojectors. Add computation of integrals of 2D monomials on the boundary of the face.

Completed first implementation of ComputeFaceProjections.

577f364

The method now fills a vector containing the integral mean of projections of basis functions whose support intersects the face.

Merge remote-tracking branch 'origin/develop' into feature/borio/vem

1ce61f9

Add computation of basis functions boundary integrals.

f117582

Add computation of piNabla projection on 3D cells.

3c68411

Bug-fixes and changes to face computations.

Add computation of the integral mean of VEM basis functions.

0578899

Completed first implementation of 3D conforming VEM.

4577a34

Add computation of integral mean of derivatives of basis functions. Add computation of the stabilization matrix. Add computation of the number of quadrature points. Code cleanups.

Merge remote-tracking branch 'origin/develop' into feature/borio/vem

6a81679

Bugfixes after merge and adding test on stabilizationMatrix.

930dc41

Face rotation matrices are now computed on-the-fly using pre-computed face normals.

Merge remote-tracking branch 'origin/develop' into feature/borio/vem

a650f07

andrea-borio marked this pull request as ready for review December 10, 2020 11:30

andrea-borio requested review from castelletto1, francoishamon, joshua-white and rrsettgast December 10, 2020 11:31

Fix bug when compiling Release.

ca1b984

andrea-borio added 2 commits December 11, 2020 10:18

Attempting to fix a bug with CUDA compilation.

0160297

Bugfix in testVirtualElementBase.

44ac67e

andrea-borio added 3 commits December 11, 2020 16:37

Avoiding arraySlices in testVirtualElementBase.cpp

e845454

Using arraySlices caused build errors with CUDA.

Merge remote-tracking branch 'origin/develop' into feature/borio/vem

b0e5a03

Applying uncrustify_style.

559a6fb

andrea-borio added the flag: ready for review label Dec 14, 2020

francoishamon approved these changes Jan 11, 2021

View reviewed changes

castelletto1 approved these changes Jan 15, 2021

View reviewed changes

andrea-borio added 4 commits January 18, 2021 21:51

Changed name of parameters of ConformingVirtualElementOrder1

5123fc3

Code cleanups following PR suggestions.

87406c9

Merge remote-tracking branch 'origin/develop' into feature/borio/vem

f4182de

Remove constructor of ConformingVirtualElementOrder1 and apply uncrus…

19575ff

…tify_style.

andrea-borio added the ci: run CUDA builds Allows to triggers (costly) CUDA jobs label Jan 19, 2021

Merge remote-tracking branch 'origin/develop' into feature/borio/vem

4ce7730

castelletto1 removed the flag: ready for review label Jan 19, 2021

castelletto1 merged commit 99613c1 into develop Jan 19, 2021

castelletto1 deleted the feature/borio/vem branch January 19, 2021 21:39

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Virtual Element Method #1220

Virtual Element Method #1220

andrea-borio commented Oct 27, 2020

andrea-borio commented Dec 10, 2020

francoishamon left a comment

francoishamon Jan 11, 2021

andrea-borio Jan 12, 2021

francoishamon Jan 11, 2021

francoishamon Jan 11, 2021

francoishamon Jan 11, 2021

francoishamon Jan 11, 2021

francoishamon Jan 11, 2021

castelletto1 Jan 15, 2021

francoishamon Jan 11, 2021

castelletto1 left a comment •

edited

Loading

castelletto1 Jan 15, 2021 •

edited

Loading

andrea-borio Jan 18, 2021

castelletto1 Jan 19, 2021

castelletto1 Jan 15, 2021

TotoGaz commented Jan 16, 2021

andrea-borio commented Jan 19, 2021 •

edited

Loading

	array2d< real64 > faceRotationMatrix( 3, 3 );
	real64 faceRotationMatrix[ 3 ][ 3 ] = {{ 0.0 }};

		array1d< real64 > monomInternalIntegrals( 3 );
		LvArray::tensorOps::fill< 3 >( monomInternalIntegrals, 0.0 );

	array1d< real64 > monomBoundaryIntegrals( 4 );
	real64 monomBoundaryIntegrals[ 4 ] = { 0.0 };

Virtual Element Method #1220

Virtual Element Method #1220

Conversation

andrea-borio commented Oct 27, 2020

andrea-borio commented Dec 10, 2020

francoishamon left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

castelletto1 left a comment • edited Loading

Choose a reason for hiding this comment

castelletto1 Jan 15, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

TotoGaz commented Jan 16, 2021

andrea-borio commented Jan 19, 2021 • edited Loading

castelletto1 left a comment •

edited

Loading

castelletto1 Jan 15, 2021 •

edited

Loading

andrea-borio commented Jan 19, 2021 •

edited

Loading