Skip to content

Commit 73e3f01

Browse files
mmalinenmmalinen
committed
A squeeze of the code
1 parent 77e1f5d commit 73e3f01

File tree

1 file changed

+33
-176
lines changed

1 file changed

+33
-176
lines changed

‎fem/src/ElementDescription.F90‎

Lines changed: 33 additions & 176 deletions
Original file line numberDiff line numberDiff line change
@@ -8701,185 +8701,39 @@ FUNCTION EdgeElementInfo( Element, Nodes, u, v, w, F, G, detF, &
87018701
END IF
87028702

87038703
IF (SecondOrder) THEN
8704-
!-------------------------------------------------
8705-
! Two basis functions defined on the face 213:
8706-
!-------------------------------------------------
8707-
TriangleFaceMap(:) = (/ 2,1,3 /)
8708-
FaceIndices(1:3) = GIndexes(TriangleFaceMap(1:3))
8709-
CALL TriangleFaceDofsOrdering(I1,I2,D1,D2,FaceIndices)
87108704

8711-
WorkBasis(1,1) = ((4.0d0*v - Sqrt(2.0d0)*w)*&
8712-
(-6.0d0 + 2.0d0*Sqrt(3.0d0)*v + Sqrt(6.0d0)*w))/(48.0d0*Sqrt(3.0d0))
8713-
WorkBasis(1,2) = -(u*(4.0d0*v - Sqrt(2.0d0)*w))/24.0d0
8714-
WorkBasis(1,3) = (u*(-2.0d0*Sqrt(2.0d0)*v + w))/24.0d0
8715-
WorkCurlBasis(1,1) = -u/(4.0d0*Sqrt(2.0d0))
8716-
WorkCurlBasis(1,2) = (Sqrt(6.0d0) + 3.0d0*Sqrt(2.0d0)*v - 3.0d0*w)/24.0d0
8717-
WorkCurlBasis(1,3) = (Sqrt(3.0d0) - 3.0d0*v)/6.0d0
8718-
8719-
WorkBasis(2,1) = ((4.0d0*v - Sqrt(2.0d0)*w)*(-6.0d0 + 6.0d0*u + &
8720-
2.0d0*Sqrt(3.0d0)*v + Sqrt(6.0d0)*w))/(96.0d0*Sqrt(3.0d0))
8721-
WorkBasis(2,2) = -((4.0d0*Sqrt(3.0d0) + 4.0d0*Sqrt(3.0d0)*u - 3.0d0*Sqrt(2.0d0)*w)*&
8722-
(-6.0d0 + 6.0d0*u + 2.0d0*Sqrt(3.0d0)*v + Sqrt(6.0d0)*w))/288.0d0
8723-
WorkBasis(2,3) = ((Sqrt(3.0d0) + Sqrt(3.0d0)*u - 3.0d0*v)*&
8724-
(-6.0d0 + 6.0d0*u + 2.0d0*Sqrt(3.0d0)*v + Sqrt(6.0d0)*w))/(144.0d0*Sqrt(2.0d0))
8725-
WorkCurlBasis(2,1) = -(-6.0d0 + 2.0d0*u + 2.0d0*Sqrt(3.0d0)*v + &
8726-
Sqrt(6.0d0)*w)/(16.0d0*Sqrt(2.0d0))
8727-
WorkCurlBasis(2,2) = (2.0d0*Sqrt(3.0d0) - 6.0d0*Sqrt(3.0d0)*u + &
8728-
6.0d0*v - 3.0d0*Sqrt(2.0d0)*w)/(48.0d0*Sqrt(2.0d0))
8729-
WorkCurlBasis(2,3) = (Sqrt(3.0d0) - 3.0d0*Sqrt(3.0d0)*u - 3.0d0*v)/12.0d0
8730-
8731-
WorkBasis(3,1) = -((4.0d0*v - Sqrt(2.0d0)*w)*(-6.0d0 - 6.0d0*u + &
8732-
2.0d0*Sqrt(3.0d0)*v + Sqrt(6.0d0)*w))/(96.0d0*Sqrt(3.0d0))
8733-
WorkBasis(3,2) = ((-4.0d0*Sqrt(3.0d0) + 4.0d0*Sqrt(3.0d0)*u + 3.0d0*Sqrt(2.0d0)*w)* &
8734-
(-6.0d0 - 6.0d0*u + 2.0d0*Sqrt(3.0d0)*v + Sqrt(6.0d0)*w))/288.0d0
8735-
WorkBasis(3,3) = -((-Sqrt(3.0d0) + Sqrt(3.0d0)*u + 3.0d0*v)* &
8736-
(-6.0d0 - 6.0d0*u + 2.0d0*Sqrt(3.0d0)*v + Sqrt(6.0d0)*w))/(144.0d0*Sqrt(2.0d0))
8737-
WorkCurlBasis(3,1) = -(-6.0d0 - 2.0d0*u + 2.0d0*Sqrt(3.0d0)*v + &
8738-
Sqrt(6.0d0)*w)/(16.0d0*Sqrt(2.0d0))
8739-
WorkCurlBasis(3,2) = (-2.0d0*Sqrt(3.0d0) - 6.0d0*Sqrt(3.0d0)*u - 6.0d0*v + &
8740-
3.0d0*Sqrt(2.0d0)*w)/(48.0d0*Sqrt(2.0d0))
8741-
WorkCurlBasis(3,3) = (-Sqrt(3.0d0) - 3.0d0*Sqrt(3.0d0)*u + 3.0d0*v)/12.0d0
8705+
DO k=1,4
8706+
!
8707+
! Two additional basis functions on each face
8708+
!
8709+
SELECT CASE(k)
8710+
CASE(1)
8711+
TriangleFaceMap(:) = (/ 2,1,3 /)
8712+
CASE(2)
8713+
TriangleFaceMap(:) = (/ 1,2,4 /)
8714+
CASE(3)
8715+
TriangleFaceMap(:) = (/ 2,3,4 /)
8716+
CASE(4)
8717+
TriangleFaceMap(:) = (/ 3,1,4 /)
8718+
END SELECT
87428719

8743-
IF (RedefineFaceBasis) THEN
8744-
EdgeBasis(13,:) = 0.5d0 * D1 * WorkBasis(I1,:) + 0.5d0 * D2 * WorkBasis(I2,:)
8745-
CurlBasis(13,:) = 0.5d0 * D1 * WorkCurlBasis(I1,:) + 0.5d0 * D2 * WorkCurlBasis(I2,:)
8746-
EdgeBasis(14,:) = 0.5d0 * D2 * WorkBasis(I2,:) - 0.5d0 * D1 * WorkBasis(I1,:)
8747-
CurlBasis(14,:) = 0.5d0 * D2 * WorkCurlBasis(I2,:) - 0.5d0 * D1 * WorkCurlBasis(I1,:)
8748-
ELSE
8749-
EdgeBasis(13,:) = D1 * WorkBasis(I1,:)
8750-
CurlBasis(13,:) = D1 * WorkCurlBasis(I1,:)
8751-
EdgeBasis(14,:) = D2 * WorkBasis(I2,:)
8752-
CurlBasis(14,:) = D2 * WorkCurlBasis(I2,:)
8753-
END IF
8754-
8755-
!-------------------------------------------------
8756-
! Two basis functions defined on the face 124:
8757-
!-------------------------------------------------
8758-
TriangleFaceMap(:) = (/ 1,2,4 /)
8759-
FaceIndices(1:3) = GIndexes(TriangleFaceMap(1:3))
8760-
CALL TriangleFaceDofsOrdering(I1,I2,D1,D2,FaceIndices)
8720+
FaceIndices(1:3) = GIndexes(TriangleFaceMap(1:3))
8721+
CALL TriangleFaceDofsOrdering(I1,I2,D1,D2,FaceIndices)
87618722

8762-
WorkBasis(1,1) = -(w*(-6.0d0 + 2.0d0*Sqrt(3.0d0)*v + Sqrt(6.0d0)*w))/(8.0d0*Sqrt(6.0d0))
8763-
WorkBasis(1,2) = (u*w)/(4.0d0*Sqrt(2.0d0))
8764-
WorkBasis(1,3) = (u*w)/8.0d0
8765-
WorkCurlBasis(1,1) = -u/(4.0d0*Sqrt(2.0d0))
8766-
WorkCurlBasis(1,2) = (Sqrt(6.0d0) - Sqrt(2.0d0)*v - 3.0d0*w)/8.0d0
8767-
WorkCurlBasis(1,3) = w/(2.0d0*Sqrt(2.0d0))
8768-
8769-
WorkBasis(2,1) = -(w*(-6.0d0 - 6.0d0*u + 2.0d0*Sqrt(3.0d0)*v + &
8770-
Sqrt(6.0d0)*w))/(16.0d0*Sqrt(6.0d0))
8771-
WorkBasis(2,2) = (w*(1.0d0 + u - v/Sqrt(3.0d0) - w/Sqrt(6.0d0)))/(8.0d0*Sqrt(2.0d0))
8772-
WorkBasis(2,3) = ((-Sqrt(3.0d0) + Sqrt(3.0d0)*u + v)* &
8773-
(-6.0d0 - 6.0d0*u + 2.0d0*Sqrt(3.0d0)*v + Sqrt(6.0d0)*w))/(48.0d0*Sqrt(2.0d0))
8774-
WorkCurlBasis(2,1) = (-3.0d0*Sqrt(2.0d0) - Sqrt(2.0d0)*u + Sqrt(6.0d0)*v + Sqrt(3.0d0)*w)/16.0d0
8775-
WorkCurlBasis(2,2) = (Sqrt(6.0d0) + 3.0d0*Sqrt(6.0d0)*u - Sqrt(2.0d0)*v - 3.0d0*w)/16.0d0
8776-
WorkCurlBasis(2,3) = w/(4.0d0*Sqrt(2.0d0))
8777-
8778-
WorkBasis(3,1) = (w*(-6.0d0 + 6.0d0*u + 2.0d0*Sqrt(3.0d0)*v + Sqrt(6.0d0)*w))/(16.0d0*Sqrt(6.0d0))
8779-
WorkBasis(3,2) = -(w*(-6.0d0 + 6.0d0*u + 2.0d0*Sqrt(3.0d0)*v + Sqrt(6.0d0)*w))/(48.0d0*Sqrt(2.0d0))
8780-
WorkBasis(3,3) = -((Sqrt(6.0d0) + Sqrt(6.0d0)*u - Sqrt(2.0d0)*v)*&
8781-
(-6.0d0 + 6.0d0*u + 2.0d0*Sqrt(3.0d0)*v + Sqrt(6.0d0)*w))/96.0d0
8782-
WorkCurlBasis(3,1) = (-3.0d0*Sqrt(2.0d0) + Sqrt(2.0d0)*u + Sqrt(6.0d0)*v + Sqrt(3.0d0)*w)/16.0d0
8783-
WorkCurlBasis(3,2) = (-Sqrt(6.0d0) + 3.0d0*Sqrt(6.0d0)*u + Sqrt(2.0d0)*v + 3.0d0*w)/16.0d0
8784-
WorkCurlBasis(3,3) = -w/(4.0d0*Sqrt(2.0d0))
8785-
8786-
IF (RedefineFaceBasis) THEN
8787-
EdgeBasis(15,:) = 0.5d0 * D1 * WorkBasis(I1,:) + 0.5d0 * D2 * WorkBasis(I2,:)
8788-
CurlBasis(15,:) = 0.5d0 * D1 * WorkCurlBasis(I1,:) + 0.5d0 * D2 * WorkCurlBasis(I2,:)
8789-
EdgeBasis(16,:) = 0.5d0 * D2 * WorkBasis(I2,:) - 0.5d0 * D1 * WorkBasis(I1,:)
8790-
CurlBasis(16,:) = 0.5d0 * D2 * WorkCurlBasis(I2,:) - 0.5d0 * D1 * WorkCurlBasis(I1,:)
8791-
ELSE
8792-
EdgeBasis(15,:) = D1 * WorkBasis(I1,:)
8793-
CurlBasis(15,:) = D1 * WorkCurlBasis(I1,:)
8794-
EdgeBasis(16,:) = D2 * WorkBasis(I2,:)
8795-
CurlBasis(16,:) = D2 * WorkCurlBasis(I2,:)
8796-
END IF
8723+
CALL WeightedWhitneyForms(WorkBasis(1:3,1:3), WorkCurlBasis(1:3,1:3), k, u, v, w)
87978724

8798-
!-------------------------------------------------
8799-
! Two basis functions defined on the face 234:
8800-
!-------------------------------------------------
8801-
TriangleFaceMap(:) = (/ 2,3,4 /)
8802-
FaceIndices(1:3) = GIndexes(TriangleFaceMap(1:3))
8803-
CALL TriangleFaceDofsOrdering(I1,I2,D1,D2,FaceIndices)
8804-
8805-
WorkBasis(1,1) = (w*(-2.0d0*Sqrt(2.0d0)*v + w))/16.0d0
8806-
WorkBasis(1,2) = (w*(4.0d0*Sqrt(3.0d0) + 4.0d0*Sqrt(3.0d0)*u - &
8807-
3.0d0*Sqrt(2.0d0)*w))/(16.0d0*Sqrt(6.0d0))
8808-
WorkBasis(1,3) = -((1.0d0 + u - Sqrt(3.0d0)*v)*w)/16.0d0
8809-
WorkCurlBasis(1,1) = (-2.0d0*Sqrt(2.0d0) - 2.0d0*Sqrt(2.0d0)*u + 3.0d0*Sqrt(3.0d0)*w)/16.0d0
8810-
WorkCurlBasis(1,2) = (-2.0d0*Sqrt(2.0d0)*v + 3.0d0*w)/16.0d0
8811-
WorkCurlBasis(1,3) = w/(2.0d0*Sqrt(2.0d0))
8812-
8813-
WorkBasis(2,1) = (w*(-2.0d0*Sqrt(2.0d0)*v + w))/16.0d0
8814-
WorkBasis(2,2) = -(w*(-4.0d0*v + Sqrt(2.0d0)*w))/(16.0d0*Sqrt(6.0d0))
8815-
WorkBasis(2,3) = -((Sqrt(6.0d0) + Sqrt(6.0d0)*u - Sqrt(2.0d0)*v)*&
8816-
(-4.0d0*v + Sqrt(2.0d0)*w))/(32.0d0*Sqrt(3.0d0))
8817-
WorkCurlBasis(2,1) = (2.0d0*Sqrt(2.0d0) + 2.0d0*Sqrt(2.0d0)*u - &
8818-
2.0d0*Sqrt(6.0d0)*v + Sqrt(3.0d0)*w)/16.0d0
8819-
WorkCurlBasis(2,2) = (-4.0d0*Sqrt(2.0d0)*v + 3.0d0*w)/16.0d0
8820-
WorkCurlBasis(2,3) = w/(4.0d0*Sqrt(2.0d0))
8821-
8822-
WorkBasis(3,1) = 0.0d0
8823-
WorkBasis(3,2) = (w*(-6.0d0 - 6.0d0*u + 2.0d0*Sqrt(3.0d0)*v + Sqrt(6.0d0)*w))/(24.0d0*Sqrt(2.0d0))
8824-
WorkBasis(3,3) = -(v*(-6.0d0 - 6.0d0*u + 2.0d0*Sqrt(3.0d0)*v + Sqrt(6.0d0)*w))/(24.0d0*Sqrt(2.0d0))
8825-
WorkCurlBasis(3,1) = (2.0d0*Sqrt(2.0d0) + 2.0d0*Sqrt(2.0d0)*u - Sqrt(6.0d0)*v - Sqrt(3.0d0)*w)/8.0d0
8826-
WorkCurlBasis(3,2) = -v/(4.0d0*Sqrt(2.0d0))
8827-
WorkCurlBasis(3,3) = -w/(4.0d0*Sqrt(2.0d0))
8828-
8829-
IF (RedefineFaceBasis) THEN
8830-
EdgeBasis(17,:) = 0.5d0 * D1 * WorkBasis(I1,:) + 0.5d0 * D2 * WorkBasis(I2,:)
8831-
CurlBasis(17,:) = 0.5d0 * D1 * WorkCurlBasis(I1,:) + 0.5d0 * D2 * WorkCurlBasis(I2,:)
8832-
EdgeBasis(18,:) = 0.5d0 * D2 * WorkBasis(I2,:) - 0.5d0 * D1 * WorkBasis(I1,:)
8833-
CurlBasis(18,:) = 0.5d0 * D2 * WorkCurlBasis(I2,:) - 0.5d0 * D1 * WorkCurlBasis(I1,:)
8834-
ELSE
8835-
EdgeBasis(17,:) = D1 * WorkBasis(I1,:)
8836-
CurlBasis(17,:) = D1 * WorkCurlBasis(I1,:)
8837-
EdgeBasis(18,:) = D2 * WorkBasis(I2,:)
8838-
CurlBasis(18,:) = D2 * WorkCurlBasis(I2,:)
8839-
END IF
8840-
8841-
!-------------------------------------------------
8842-
! Two basis functions defined on the face 314:
8843-
!-------------------------------------------------
8844-
TriangleFaceMap(:) = (/ 3,1,4 /)
8845-
FaceIndices(1:3) = GIndexes(TriangleFaceMap(1:3))
8846-
CALL TriangleFaceDofsOrdering(I1,I2,D1,D2,FaceIndices)
8847-
8848-
WorkBasis(1,1) = (w*(-2.0d0*Sqrt(2.0d0)*v + w))/16.0d0
8849-
WorkBasis(1,2) = (w*(-4.0d0*Sqrt(3.0d0) + 4.0d0*Sqrt(3.0d0)*u + &
8850-
3.0d0*Sqrt(2.0d0)*w))/(16.0d0*Sqrt(6.0d0))
8851-
WorkBasis(1,3) = -((-1.0d0 + u + Sqrt(3.0d0)*v)*w)/16.0d0
8852-
WorkCurlBasis(1,1) = (2.0d0*Sqrt(2.0d0) - 2.0d0*Sqrt(2.0d0)*u - 3.0d0*Sqrt(3.0d0)*w)/16.0d0
8853-
WorkCurlBasis(1,2) = (-2.0d0*Sqrt(2.0d0)*v + 3.0d0*w)/16.0d0
8854-
WorkCurlBasis(1,3) = w/(2.0d0*Sqrt(2.0d0))
8855-
8856-
WorkBasis(2,1) = 0.0d0
8857-
WorkBasis(2,2) = (w*(-6.0d0 + 6.0d0*u + 2.0d0*Sqrt(3.0d0)*v + Sqrt(6.0d0)*w))/(24.0d0*Sqrt(2.0d0))
8858-
WorkBasis(2,3) = -(v*(-6.0d0 + 6.0d0*u + 2.0d0*Sqrt(3.0d0)*v + Sqrt(6.0d0)*w))/(24.0d0*Sqrt(2.0d0))
8859-
WorkCurlBasis(2,1) = (2.0d0*Sqrt(2.0d0) - 2.0d0*Sqrt(2.0d0)*u - Sqrt(6.0d0)*v - Sqrt(3.0d0)*w)/8.0d0
8860-
WorkCurlBasis(2,2) = v/(4.0d0*Sqrt(2.0d0))
8861-
WorkCurlBasis(2,3) = w/(4.0d0*Sqrt(2.0d0))
8862-
8863-
WorkBasis(3,1) = ((2.0d0*Sqrt(2.0d0)*v - w)*w)/16.0d0
8864-
WorkBasis(3,2) = -(w*(-4.0d0*v + Sqrt(2.0d0)*w))/(16.0d0*Sqrt(6.0d0))
8865-
WorkBasis(3,3) = ((-Sqrt(3.0d0) + Sqrt(3.0d0)*u + v)*&
8866-
(-4.0d0*v + Sqrt(2.0d0)*w))/(16.0d0*Sqrt(6.0d0))
8867-
WorkCurlBasis(3,1) = (2.0d0*Sqrt(2.0d0) - 2.0d0*Sqrt(2.0d0)*u - &
8868-
2.0d0*Sqrt(6.0d0)*v + Sqrt(3.0d0)*w)/16.0d0
8869-
WorkCurlBasis(3,2) = (4.0d0*Sqrt(2.0d0)*v - 3.0d0*w)/16.0d0
8870-
WorkCurlBasis(3,3) = -w/(4.0d0*Sqrt(2.0d0))
8871-
8872-
IF (RedefineFaceBasis) THEN
8873-
EdgeBasis(19,:) = 0.5d0 * D1 * WorkBasis(I1,:) + 0.5d0 * D2 * WorkBasis(I2,:)
8874-
CurlBasis(19,:) = 0.5d0 * D1 * WorkCurlBasis(I1,:) + 0.5d0 * D2 * WorkCurlBasis(I2,:)
8875-
EdgeBasis(20,:) = 0.5d0 * D2 * WorkBasis(I2,:) - 0.5d0 * D1 * WorkBasis(I1,:)
8876-
CurlBasis(20,:) = 0.5d0 * D2 * WorkCurlBasis(I2,:) - 0.5d0 * D1 * WorkCurlBasis(I1,:)
8877-
ELSE
8878-
EdgeBasis(19,:) = D1 * WorkBasis(I1,:)
8879-
CurlBasis(19,:) = D1 * WorkCurlBasis(I1,:)
8880-
EdgeBasis(20,:) = D2 * WorkBasis(I2,:)
8881-
CurlBasis(20,:) = D2 * WorkCurlBasis(I2,:)
8882-
END IF
8725+
IF (RedefineFaceBasis) THEN
8726+
EdgeBasis(12+2*(k-1)+1,:) = 0.5d0 * D1 * WorkBasis(I1,:) + 0.5d0 * D2 * WorkBasis(I2,:)
8727+
CurlBasis(12+2*(k-1)+1,:) = 0.5d0 * D1 * WorkCurlBasis(I1,:) + 0.5d0 * D2 * WorkCurlBasis(I2,:)
8728+
EdgeBasis(12+2*k,:) = 0.5d0 * D2 * WorkBasis(I2,:) - 0.5d0 * D1 * WorkBasis(I1,:)
8729+
CurlBasis(12+2*k,:) = 0.5d0 * D2 * WorkCurlBasis(I2,:) - 0.5d0 * D1 * WorkCurlBasis(I1,:)
8730+
ELSE
8731+
EdgeBasis(12+2*(k-1)+1,:) = D1 * WorkBasis(I1,:)
8732+
CurlBasis(12+2*(k-1)+1,:) = D1 * WorkCurlBasis(I1,:)
8733+
EdgeBasis(12+2*k,:) = D2 * WorkBasis(I2,:)
8734+
CurlBasis(12+2*k,:) = D2 * WorkCurlBasis(I2,:)
8735+
END IF
8736+
END DO
88838737

88848738
! Finally, scale to reduce ill-conditioning:
88858739
IF (ScaleFaceBasis) THEN
@@ -9640,6 +9494,9 @@ FUNCTION EdgeElementInfo( Element, Nodes, u, v, w, F, G, detF, &
96409494
EdgeMap => GetEdgeMap(7)
96419495

96429496
IF (SecondOrder) THEN
9497+
! IF (Simplicial) THEN
9498+
9499+
! ELSE
96439500
!---------------------------------------------------------------
96449501
! The second-order element from the Nedelec's first family
96459502
! (note that the lowest-order prism element is from a different
@@ -10015,7 +9872,7 @@ FUNCTION EdgeElementInfo( Element, Nodes, u, v, w, F, G, detF, &
100159872
EdgeBasis(35:36,1:2) = sqrt(150.0d0) * EdgeBasis(35:36,1:2)
100169873
CurlBasis(35:36,1:3) = sqrt(150.0d0) * CurlBasis(35:36,1:3)
100179874
END IF
10018-
9875+
! END IF
100199876
ELSE
100209877
!--------------------------------------------------------------
100219878
! The lowest-order element from the optimal family. The optimal

0 commit comments

Comments
 (0)