# Spatio-Temporal Video Representation Learning for
AI Based Video Playback Style Prediction

Rishubh Parihar
Indian Institute of Science,
Bangalore, India.
<EMAIL_ADDRESS>Equal Contribution Gaurav Ramola 11footnotemark: 1
Samsung India Research Institute,
Bangalore, India.
<EMAIL_ADDRESS>Ranajit Saha
Microsoft Corporation,
Hyderabad, India.
<EMAIL_ADDRESS>Ravi Kini
Samsung India Research Institute,
Bangalore, India.
<EMAIL_ADDRESS>Aniket Rege
Univ. of Washinton,
Seattle, USA.
<EMAIL_ADDRESS>Sudha Velusamy
Samsung India Research Institute,
Bangalore, India.
<EMAIL_ADDRESS>
###### Abstract

Ever-increasing smartphone-generated video content demands intelligent
techniques to edit and enhance videos on power-constrained devices. Most of
the best performing algorithms for video understanding tasks like action
recognition, localization, etc., rely heavily on rich spatio-temporal
representations to make accurate predictions. For effective learning of the
spatio-temporal representation, it is crucial to understand the underlying
object motion patterns present in the video. In this paper, we propose a novel
approach for understanding object motions via motion type classification. The
proposed motion type classifier predicts a motion type for the video based on
the trajectories of the objects present. Our classifier assigns a motion type
for the given video from the following five primitive motion classes: linear,
projectile, oscillatory, local and random. We demonstrate that the
representations learned from the motion type classification generalizes well
for the challenging downstream task of video retrieval. Further, we proposed a
recommendation system for video playback style based on the motion type
classifier predictions.

##  1 Introduction

An increasing volume of smart-phones with high-quality cameras in recent years
has led to a meteoric rise in the amount of video content captured and shared
on social media platforms such as Tiktok, YouTube, Facebook, Instagram,
SnapChat, ShareChat etc. This trend has fostered the need for automated video
analysis tools that can aid the user to edit videos with ease on mobile
devices, on-the-fly.

Figure 1: Visualizing an example motion trajectory of a ball

Videos contain rich information embedded in both spatial and temporal
dimensions, which together capture the overall dynamics of the scene. Learning
meaningful spatio-temporal representation is at the core of most video
analysis tasks like video retrieval, action recognition, temporal and spatial
action localization, object motion analysis, video captioning, and modelling
of human-object interactions. There is a fundamental need for methods to learn
generalized spatio-temporal representations that can work effectively for
multiple downstream tasks. One of the popular approaches is to train a model
for video action recognition and obtain the implicitly learned video
representation [3] [31]. Recently, many self-supervised methods have been
proposed, where a deep network is trained for an auxiliary pre-text task to
learn rich spatio-temporal representations.

Object motion understanding is crucial to learn rich spatio-temporal
representations as it provides insights about the natural motion pattern of
objects in the world and how they interact with other objects in the scene
[38]. For instance, consider the example of a video where a person is shooting
a ball towards the goalpost as shown in Fig. 1. Analysing the motion of the
ball during this action will provide insight about the most likely motion of
the soccer ball: just after kicking, the ball will follow a projectile motion
in the air, and after dropping on floor the ball will bounce a few times. This
motion pattern of a relatively common occurrences in everyday life is
extremely complex to model in a mathematical or mechanical sense as it
comprises, for instance in the above example, movement of the player’s body
and real world forces (friction, air drag) at play.

In this work, we present a method of analysing the underlying directional
information of object motions that occur in real-world human actions like
kicking, walking, jumping, clapping, etc., by estimating the object motion
type in a video. As it is difficult to jointly model motions of all the
objects in the scene, we focus only on the dominant motion in the video. To
this end, we have formulated a classification problem, to classify the
directional motion pattern into one of the defined classes. Based on our
internal study on action classes present in popular video dataset HMDB$51$
[21], we have defined five primitive motion type classes: _linear, projectile,
oscillatory, local and random._ According to us, most of the real world human
actions can be assigned to one of the above defined motion types. For
instance: _walking, running, bike-riding_ have a linear motion type as the
dominant motion, _kicking, cartwheel_ makes projectile motion, and _talking,
chewing, smoking_ have a local motion type. All the motion patterns having
periodic motion are considered under oscillatory class, for example, _pushup
and exercise._ The actions which do not lie into any of these categories were
assigned the class random. To our knowledge, there is no open-source video
dataset currently available with motion type labels for videos. To this end,
we have added motion type annotations to the HMD51 [21] dataset for training
the motion classifier. The motivation of this work is to address the
following:_1) Is it possible for a neural network model to perform well on the
task of motion type classification? 2) What internal feature representations
does the model learn in this process? 3) Are these learned features generalize
well on other downstream video understanding tasks?_ We have tried to answer
these questions throughout this paper by training a CNN model for motion type
classification and analyzing its learned features through general video
analysis tasks like video retrieval.

We also demonstrate an exciting use-case of the above-presented motion type
classification method: video playback style recommendation, which boosts the
overall aesthetics of the videos. A few common playback styles include:
Reverse (temporally reversing the video), Loop (repeating the video in a
loop), Boomerang (playing a concatenated video of normal and reverse). Finding
a suitable playback style is often a time-consuming process where a user
manually applies each available playback style. This created a space to
engineer automated tools for this problem. Our proposed solution tries to
automate this process of playback style selection. More details for the design
of this recommendation algorithm are presented in Sec. 3.2.

Lastly, we show that through the proposed motion type classification, we are
able to learn rich spatio-temporal representations that generalize well for
other video analysis tasks such as video retrieval. In a subjective evaluation
of the learned representations for video retrieval, we achieved promising
results on the HMDB$51$ dataset. Furthermore, we made specific design choices
to make the network efficient for mobile deployment. Our model for motion
classification has inference time of $200ms$ for a $10$ second video clip on a
Samsung S20 phone.

We summarize our major contributions as follows:

  1. 1.

A neural network for understanding object motion in videos by classifying
object motion type into one of the five primitive motion classes: _linear,
projectile, oscillatory, local and random._

  2. 2.

A light-weight network for video representation learning that is suitable for
real-time execution on mobile devices.

  3. 3.

A recommender system to predict suitable video playback style for videos by
analysing predicted object motion patterns

##  2 Related Works

Figure 2: Overall network architecture for motion type classification. Given
an input video, we divide it temporally into three segments and extract one
central frame from each segment. These three frames are fed to the Feature
Extractor Network, and the extracted features are then averaged to obtain a
$1280$-dimensional (1280D) feature vector, which is used for motion type
classification

Video action recognition has been studied extensively by computer vision
community. The success of video action recognition majorly depends on crafting
the spatio-temporal features in the video representations. Traditionally,
video features are extracted from optical-flow based motion information in the
videos, e.g. Motion Boundary Histograms (MBH) [5] and trajectories [32],
Histograms Of Flow (HOF) [22] or spatiotemporal oriented filtering such as
HOG3D [20], Cuboids [8] and Spatiotemporal Oriented Energies (SOEs) [11, 7].
The resounding success of Convolutional Neural Networks (CNNs) for image
processing applications has caused its extension to video processing problems
as well. Just like the spatial features, deep CNNs are also capable to extract
accurate temporal information as well e.g. FlowNets [9, 17]. Both the temporal
and spatial information are important in various video recognition tasks.
Simonyan and Zisserman [29] has proposed a two-stream CNN architecture to
incorporate both spatial and temporal features of the videos. The spatial
features are captured by passing the RGB frames of the videos and the temporal
features are captured by extracting the flow frames. Several other works [12,
13] have explored the different effective fusion options of two streams - flow
and RGB streams. The major bottleneck in two-stream networks as well as
optical flow based methods is the optical flow extraction step as it consumes
a lot of time and hence the inference time increases.

DMC-Net [28] approximates the flow using a reconstruction loss and an
adversarial loss jointly for the task of action classification. This model is
two folds faster than the state-of-the-art methods and achieves accuracy close
to the methods using optical flow information. The study of Tran _et al_. [30]
shows the effectiveness of using $3D$-CNNs instead of $2D$-CNNs to model both
spatial and temporal features together in a single branch. Although $3D$-CNNs
produce promising results, it is much more expensive than $2D$-CNNs.
Experiments by Xie _et al_. [36] showed that we can trade-off accuracy and
speed by replacing some $3D$ conv layers by $2D$ convolutions. Having $3D$
conv layers at the higher layers and $2D$ conv layers at the lower part of the
network is faster and this configuration surprisingly has higher accuracy.
They also propose separable $3D$-CNN (S3D) configuration which separates
spatial and temporal $3D$ convolutions. MARS [4] introduces the learning
approaches to train $3D$-CNN operating on RGB frames which mimics the motion
stream. It eradicates the need of flow extraction during the inference time.
Frame sampling from videos is also an important part in video processing.
Temporal Segment Network (TSN) [34] works on sparse temporal snippets. The
videos are split into k chunks and a small snippet in chosen from each of the
chunk. The chunks are processed individually and at the end the decisions are
aggregated as per the consensus function to come to the final conclusion. TSN
gives promising result for action recognition task. Lin _et al_. [23] proposes
a generic module called Temporal Shift module (TSM). It is a ”plug and play”
module in a network designed for video understanding task. TSM has high
efficiency and high performance. It maintains the complexity of $2D$-CNN and
performance of the $3D$-CNN. TSM facilitates the information exchange by
shifting a part of the channels along temporal dimension.

Object motion pattern understanding is crucial for learning strong spatio-
temporal features for downstream video analysis tasks [38]. There are
approaches which try to capture the object motions in the videos via learning
flow features from the videos [10, 25]. These methods predict pixel-level
feature maps for every time frame in the video, which essentially captures
only local motion patterns.

Most of the methods discussed above are based on the supervised learning
technique. But due to the scarcity of publicly available labeled dataset, it
is difficult to train deep networks with supervised learning. Several Self-
supervised methods [2, 15] for video tasks have been studied by the computer
vision community. Qian proposed [26] self-supervised Contrastive Video
Representation Learning (CVRL) method which uses the contrastive loss to map
the video clips in the embedding space. It is desired that in the embedding
space the distance between two clips from the same video is lesser than the
clips from different videos. Jenni _et al_. [18] introduced a novel self-
supervised framework to learn video representations which are sensitive to the
changes in the motion dynamics. They have observed that the motion of objects
is essential for action recognition tasks. In the proposed work, we build on
the above intuition to show that a deep network can learn rich representations
by training for motion classification.

##  3 Methodology

Humans largely use primary motion cues like underlying object motion patterns
to understand video semantics like actions or events in a scene. To perform
well on video analysis tasks like action recognition and localization, the
motion pattern representations require a semantic understanding of both the
appearance and dynamics features of the video. We aim to learn rich spatio-
temporal video representations through classification of the motion type based
on the directional motion information present in the video. To this end, we
trained a motion type classification model that classifies a video into one of
the following five primitive classes we define: _linear, oscillatory, local,
projectile, and random_. We observed that the trajectories of most natural
object motions that we encounter in the real-world can be categorized into the
first four motion classes. As it is difficult to jointly model motions of all
the objects in the scene, we focus only on the dominant motion in the video.
For instance, actions such as _walk_ and _run_ usually follow a linear
trajectory and have a dominant linear motion. Many activities that we perform
indoors have motion in only small local regions like _eat, drink, chew, talk_.
Some of the examples of actions having dominant oscillatory motion type are
_dribble, cartwheel and sit-up_. _Catch, throw, golf_ are examples for
dominant projectile motion type. Actions which do not follow any of these
directional patterns, are considered random, for instance _dance and fight_.
Some of the common real-world actions and their corresponding motion types are
shown in Table 1. To validate the quality of our learned representations, we
used these representations for video retrieval task as explained in Sec. 4.4.
As there is no publicly available video dataset with motion type labels, we
have annotated the HMDB$51$ dataset with motion type labels to obtain mHMDB51
dataset as seen in Sec. 4.1. The core of our method is a Deep Convolutional
Neural Network (Fig. 2), which is trained in a supervised fashion on mHMDB$51$
dataset for a five class motion-type classification problem.

Table 1: Mapping of real world actions to motion type and Video Playback Style in the mHMDB$51$ dataset. Example Action | Motion Type | Playback Style
---|---|---
Walk, Run | Linear | Reverse
Dive, Throw | Projectile | Boomerang
Eat, Clap | Local | Loop
PullUp | Oscillatory | Loop
Dance, Fight | Random | Random

###  3.1 Network Architecture

Most state-of-the-art networks for video representation learning and action
recognition methods [16] [3] rely on $3D$ convolutions due to their ability to
jointly learn both spatial and temporal features. However, $3D$ convolutions
have significantly higher computational cost than $2D$ convolutions, which
make them unsuitable for mobile applications that have strict power and
latency constraints. Our network uses a backbone of only $2D$ convolutions
with added Temporal Shift Modules (TSM) [23] to facilitate an information
exchange between temporally adjacent frames. This results in a light-weight
network architecture that needs very limited memory and compute requirement.
The proposed network architecture is shown in Fig.2. Our network is inspired
by TSN [34] architecture, where a set of frames is sampled from a video and
processed independently. Finally, a consensus is taken to obtain a global
feature representation. We first divide the input video temporally into $T$
segments of equal durations and one representative central frame is sampled
from each segment. The input of our model is thus a $T*N*N$ volume, where $T$
is the number of segments from the video and $N$ is both the height and the
width of the video. The input volume is passed through a TSN-style backbone
network to obtain a $T*1280$ shape feature representation. The obtained
feature vector is then averaged over the temporal dimension to obtain a
combined 1280-dimension feature vector for the entire video. This global video
feature vector is then fed into a classifier head having two fully connected
layers with $128$ and $64$ neurons respectively, followed by a softmax layer
for classification. The working of the original TSN architecture is explained
by the equation 1. The video $V$ is divided into $K$ segments
{$S_{1}$,$S_{2}$, …, $S_{K}$} of equal duration and ($T_{1}$,$T_{2}$, …,
$T_{K}$) are the sequence of snippets where each $T_{K}$ is sampled from its
corresponding segment $S_{K}$. $\mathcal{F}$($T_{K}$; $W$) defines the output
after passing the snippet $T_{K}$ through the ConvNet with parameters $W$. The
consensus module $\mathcal{G}$ combines the extracted features of all the
snippets through $\mathcal{F}$ operation. The consensus module for our
architecture takes the average of the features. The average output of
consensus module is passed through the fully-connected layer with a softmax at
the end to get the final class label. This operation is defined by
$\mathcal{H}$ in the equation 1.

$\begin{split}TSN(T_{1},T_{2},\dots,T_{K})=\mathcal{H}(\mathcal{G}(\mathcal{F}(T_{1};W),\\\
\mathcal{F}(T_{2};W),\dots,\mathcal{F}(T_{K};W)))\end{split}$ (1)

We added TSM modules in the backbone network to help the network learn strong
temporal relations across segments via shifting of the intermediate feature
channels of one segment to neighboring segments. To further reduce the
computational complexity of our network, we have used MobileNetV2 [27] as the
backbone due to its low computational cost. Our specific design choices for
the network architecture makes it suitable for video processing on mobile
devices having low compute budget.

###  3.2 Video Playback Style Recommendation

Applying a suitable playback style to a video can enhance a video and make it
more likely to be shared. Motion patterns present in the videos play an
important role in selecting the most suited playback style for the video. For
instance, for a video having linear motion like running, applying Boomerang
type will make the video counter intuitive and hence interesting. To this end,
we have designed a system for video playback style recommendation based on
predictions from motion type classifier. We have considered three most widely
used playback styles for recommendation namely Boomerang, Loop and Reverse.
Specifically, we have introduced a mapping from motion type to a suitable
playback style for an input video based on a user survey of 14 volunteers. In
this study, we showed a few example actions for each motion type to each
volunteer and asked them to select the best-suited playback style for that
corresponding action. We aggregated the results from each volunteer and
selected the most voted playback style for each motion type for the mapping.
From the results of the study as shown in Table 1, we observe that the Reverse
effect suits linear actions, and projectile motion looks good with a Boomerang
effect. For both oscillatory and local motion, loop is the best-suited
playback style. For random motion type, we randomly apply Boomerang, Reverse
or Loop. We have performed a subjective study for evaluation of our video
playback style recommendation system which is detailed in Sec. 4.3.

##  4 Experiments

We have done multiple experiments for comprehensive evaluation of our proposed
motion type classifier model. In Section 4.2 we perform an ablation with
various pre-trained weights to examine the impact of weight initialization. To
evaluate the quality of the learnt representations through motion
classification, we have performed video retrieval as detailed in Sec. 4.4. We
have also performed a subjective study for evaluation of our video playback
recommendation system. To prepare our training data, each video was first
resized: the smaller dimension was set to $256$ pixels wide, and a random
square region was cropped of side length $d$ $\epsilon$ $(256,224,192,169)$ ,
followed by a random horizontal flip. Finally, the crop was resized to
$(256,256)$ and the pixel values were normalized to the range $(0,1)$. In the
testing phase, we resized the smaller dimension to $256$ and took a center
crop. We used $T=3$ segments in all of our experiments unless mentioned
otherwise, and sampled the temporally central frame from each segment. These
three frames are the input to the network. For training, we used an initial
learning rate of $0.001$ and a learning rate schedule to reduce the learning
rate by half after the $20^{th}$ and $40^{th}$ epoch. The network was trained
for a total of $200$ epochs. Stochastic Gradient Descent was used for
optimization with momentum value of $0.9$ and a weight decay of $5e-5$. We
have trained all our models with a single P100 GPU and each training
configuration took 4hrs to converge.

###  4.1 Dataset

For all our experiments, we use the HMDB$51$ [21] dataset. The HMDB$51$
dataset contains short videos (1-15 seconds) for $51$ human actions like
cycling, eating, running and dancing etc. We have used the split-1 set of
HMDB$51$ provided by [21] to create the train/test/validation set. There are
$3570$ videos in the train set, $1530$ videos in the test set and $1749$ in
the validation set. These videos are collected from YouTube and digitized
movies and have large variability in camera motion, view-point and
illumination. For our purpose, we annotated each of the $51$ action classes
from the HMDB$51$ dataset with one of our five defined motion types. We have
named this annotated version of the HMDB$51$ dataset the mHMDB$51$ dataset. A
subset of this mapping is shown in Table 1, while the full version can be
found in the appendix.

###  4.2 Motion Classifier

For evaluation, we have compared our model with a optical flow based baseline
model and performed an ablation study with various pre-training methods. The
results are shown in Table 2.

####  4.2.1 Baseline Classifier

To benchmark our motion type classifier, we designed a baseline classifier as
a two-layer fully connected neural network. The input to this classifier is
based on the statistics of motion magnitudes in the video. To extract the
input features, we first compute the pixel-wise average over time of the
motion boundaries for the input video and divide it into $16$ cells as in
[33]. We use the standard deviation of the magnitude of motion boundaries
within each cell to form the 16-dimensional input feature vector to the motion
type classifier. In the network design, there are $128$ neurons in the first
hidden layer and $5$ neurons in the second hidden layer for the network. ReLU
activation was used after the first hidden layer and a softmax activation was
applied after the second hidden layer for the final classification. Dropout
regularization was applied with a drop probability of $0.2$ and the classifier
was trained for $5$ epochs with a learning rate of $0.001$.

####  4.2.2 Model Performance Analysis

Figure 3: Playback style recommendation by our system for YouTube videos

We observed that training our classifier from scratch achieved a performance
boost of nearly $13\%$ over the baseline flow-based model, but still low as
compared to fully supervised pre-training with ImageNet [6] and Kinetics [19].
This was expected behavior, as our model was trained with only $3500$ videos
from the HMDB$51$ dataset, which is insufficient for supervised training when
compared to the millions of data points used to train existing ImageNet and
Kinetics classifiers. Thus our usage of transfer learning via initializing our
classifier with weights learned from the ImageNet classification task
increased our accuracy by a margin of around $14\%$, due to the pre-trained
understanding of important spatial features. Initializing with weights learned
for action classification on the Kinetics dataset achieved the best accuracy,
as they have a pre-trained understanding of both spatial and temporal
features, which are useful to perform motion classification. Our baseline
local-flow-based classifier expectedly performed the worst. These observations
indicate that accurately predicting object motion type requires global
semantic information contained in motion patterns. our results demonstrate
that our motion type classifier learns more than just the motion magnitude,
and has a deeper understanding of object motion patterns.

####  4.2.3 Model Complexity Analysis

We also performed an ablation study by varying the complexity of the backbone
network. Our baseline model is TSN [34] with shift modules which process
multiple segments from a video and fuse them together at a later stage to
obtain the combined feature vector. As the number of segments represent the
complexity of the model, we have trained models with 1, 2, 3 and 8 segments in
our ablation study. The overall accuracy of the model and the number of
multiply-accumulate (MAC) operations is shown in Table 3. The three-segment
model achieved the best accuracy for motion type classification. However, the
two-segment model was able to achieve comparable accuracies to the three-
segment model with just 0.82G MAC operations, making it the optimally suited
configuration for mobile deployment. The inference time for the two-segment
model on a Samsung S20 mobile device running a Qualcomm Snapdragon Adreno 650
GPU is just 200 milliseconds. The single-segment model processes only a single
frame from the complete video and therefore struggles to learn temporal
dynamics of the video. However, it was still able to achieve a reasonable
accuracy of $61.75\%$ for motion type classification, demonstrating the
importance of object appearance in determining the natural motion patterns for
an object. The eight-segment model did not perform well, due to the HMDB$51$
dataset having small action videos and thus not requiring too many frames for
effective motion pattern understanding. We believe that passing a large number
of frames for actions with short duration captures multiple motion types
present in the video at different instances and hence confuses the network
training.

Table 2: Motion Type Classifier Top-1 Accuracy. Method | Accuracy
---|---
Baseline Classifier | 25.64
$\text{Ours}_{\text{Scratch}}$ | 38.56
$\text{Ours}_{\text{ImageNet}}$ | 57.58
$\text{Ours}_{\text{Kinetics}}$ | 72.68
Figure 4: Subjective study for video playback style recommendation.

###  4.3 Video Playback Style Recommendation

Table 3: Comparison of motion classifier top-1 accuracy and MAC operations for a varying number of input segments for the network. Segments | Accuracy | MACs
---|---|---
1 | 61.76 | 0.41G
2 | 71.05 | 0.82G
3 | 72.68 | 1.23G
8 | 68.17 | 3.28G

For a subjective evaluation of our video playback style recommendations, we
conducted a user study with $10$ volunteers. We downloaded two clips for each
of the following five actions from YouTube: _cartwheel, diving, running,
clapping, and drinking_. Our network predicted the motion type of each video,
and we applied the matching playback style based on the mapping shown in Table
1. We also prepared a comparison set for the same videos with randomly applied
playback styles. We evaluated our recommended playback styles against these
randomly selected playback styles. The volunteers were asked to select the
most aesthetic and preferred result from these two sets, the results of which
are shown in Fig. 4. For the categories that have a large global motions like
cartwheel, diving, and running, our predicted playback style was ranked better
than random playback style on an average. To our surprise, while the diving
action was not present in our training set, our engine was able to recommend
the best-suited playback style for the class. This provides evidence for the
proposition that training to predict motion type captures more abstract
information than actions, and generalizes well for unseen data. On the
contrary, for local action categories such as drinking and clapping, our
method was indistinguishable to random selection as the impact of playback
style is not very evident when motion is confined to a small spatial region.

###  4.4 Video Retrieval

To further analyze the spatio-temporal features learned by our motion type
classifier, we used these features to perform video retrieval. Given a query
video, we aim to find the three most similar videos to the query video from a
database of videos. We feed all the videos from HMDB$51$ to our motion
classifier and extract the $1280$-dimensional feature vector described in Sec.
3.1 for each video. In an ideal scenario, this feature vector represents the
motion present in the video in a compressed form. We apply the k-nearest-
neighbor algorithm in the $1280$-dimensional feature vector space to find
videos having similar motion patterns as that of the query video. Some example
retrievals from HMDB$51$ are shown in Fig. 5, from which it is evident that
our learned representations capture meaningful semantic information of object
motion. In Fig. 5a) the query video was of smoking, and all retrieved results
(laugh, chew and chew) have local facial motions. In Fig. 5b) and c) the first
two results are from the same scene but at different points in time. In Fig.
5b) the third retrieved result is of a golf swing, which has similar hand
movement to that of a cartwheel. Similarly, for c) the last retrieved result
is of a person diving from a cliff, which is very similar to the query video
of a goalkeeper diving for football. For Fig. 5d) all retrieved videos have
linear motion and in Fig. 5e) all the retrieved actions for the query video of
throw follow projectile motion.

Figure 5: Video retrieval from HMDB51 dataset using learned feature vectors

##  5 Conclusion

In this work, we have examined the importance of object motion features in
video analysis. We trained a model that understands the underlying object
motion patterns and classifies the object motion into one of the five defined
directional motion classes. We have also shown the exciting use case of
playback style recommendation based on our classifier’s predicted motion type.
Finally, we have evaluated the representations learned by motion type
classifiers for video retrieval and have found that these representations
generalize well for this task. In the future, we plan to explore other
possible approaches to model object motions in the videos. We will also
evaluate the generalization ability of learned representations for more
challenging video tasks such as action localization and classification.

## References

  * [1] Nadine Behrmann, Jurgen Gall, and Mehdi Noroozi.  Unsupervised video representation learning by bidirectional feature prediction.  In Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision, pages 1670–1679, 2021.
  * [2] Sagie Benaim, Ariel Ephrat, Oran Lang, Inbar Mosseri, William T Freeman, Michael Rubinstein, Michal Irani, and Tali Dekel.  Speednet: Learning the speediness in videos.  In Proceedings of CVPR, pages 9922–9931, 2020.
  * [3] Joao Carreira and Andrew Zisserman.  Quo vadis, action recognition? a new model and the kinetics dataset.  In proceedings of CVPR, pages 6299–6308, 2017.
  * [4] Nieves Crasto, Philippe Weinzaepfel, Karteek Alahari, and Cordelia Schmid.  Mars: Motion-augmented rgb stream for action recognition.  In Proceedings of CVPR, pages 7882–7891, 2019.
  * [5] Navneet Dalal, Bill Triggs, and Cordelia Schmid.  Human detection using oriented histograms of flow and appearance.  In ECCV, pages 428–441. Springer, 2006.
  * [6] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In 2009 CVPR, pages 248–255. Ieee, 2009.
  * [7] Konstantinos G Derpanis, Mikhail Sizintsev, Kevin J Cannons, and Richard P Wildes.  Action spotting and recognition based on a spatiotemporal orientation analysis.  IEEE transactions on pattern analysis and machine intelligence, 35(3):527–540, 2012.
  * [8] Piotr Dollár, Vincent Rabaud, Garrison Cottrell, and Serge Belongie.  Behavior recognition via sparse spatio-temporal features.  In 2005 IEEE International Workshop on Visual Surveillance and Performance Evaluation of Tracking and Surveillance, pages 65–72. IEEE, 2005\.
  * [9] Alexey Dosovitskiy, Philipp Fischer, Eddy Ilg, Philip Hausser, Caner Hazirbas, Vladimir Golkov, Patrick Van Der Smagt, Daniel Cremers, and Thomas Brox.  Flownet: Learning optical flow with convolutional networks.  In Proceedings of ICCV, pages 2758–2766, 2015.
  * [10] Lijie Fan, Wenbing Huang, Chuang Gan, Stefano Ermon, Boqing Gong, and Junzhou Huang.  End-to-end learning of motion representation for video understanding.  In Proceedings of CVPR, pages 6016–6025, 2018.
  * [11] Christoph Feichtenhofer, Axel Pinz, and Richard P Wildes.  Dynamically encoded actions based on spacetime saliency.  In Proceedings of CVPR, pages 2755–2764, 2015.
  * [12] Christoph Feichtenhofer, Axel Pinz, and Richard P Wildes.  Spatiotemporal multiplier networks for video action recognition.  In Proceedings of CVPR, pages 4768–4777, 2017.
  * [13] Christoph Feichtenhofer, Axel Pinz, and Andrew Zisserman.  Convolutional two-stream network fusion for video action recognition.  In Proceedings of CVPR, pages 1933–1941, 2016.
  * [14] Basura Fernando, Hakan Bilen, Efstratios Gavves, and Stephen Gould.  Self-supervised video representation learning with odd-one-out networks.  In Proceedings of CVPR, pages 3636–3645, 2017.
  * [15] Tengda Han, Weidi Xie, and Andrew Zisserman.  Video representation learning by dense predictive coding.  In Proceedings of CVPRW, pages 0–0, 2019.
  * [16] Kensho Hara, Hirokatsu Kataoka, and Yutaka Satoh.  Can spatiotemporal 3d cnns retrace the history of 2d cnns and imagenet?  In Proceedings of the CVPR, pages 6546–6555, 2018.
  * [17] Eddy Ilg, Nikolaus Mayer, Tonmoy Saikia, Margret Keuper, Alexey Dosovitskiy, and Thomas Brox.  Flownet 2.0: Evolution of optical flow estimation with deep networks.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2462–2470, 2017.
  * [18] Simon Jenni, Givi Meishvili, and Paolo Favaro.  Video representation learning by recognizing temporal transformations.  arXiv preprint arXiv:2007.10730, 2020.
  * [19] Will Kay, Joao Carreira, Karen Simonyan, Brian Zhang, Chloe Hillier, Sudheendra Vijayanarasimhan, Fabio Viola, Tim Green, Trevor Back, Paul Natsev, et al.  The kinetics human action video dataset.  arXiv preprint arXiv:1705.06950, 2017.
  * [20] Alexander Klaser, Marcin Marszałek, and Cordelia Schmid.  A spatio-temporal descriptor based on 3d-gradients.  2008\.
  * [21] Hildegard Kuehne, Hueihan Jhuang, Estíbaliz Garrote, Tomaso Poggio, and Thomas Serre.  Hmdb: a large video database for human motion recognition.  In 2011 ICCV, pages 2556–2563. IEEE, 2011.
  * [22] Ivan Laptev, Marcin Marszalek, Cordelia Schmid, and Benjamin Rozenfeld.  Learning realistic human actions from movies.  In 2008 CVPR, pages 1–8. IEEE, 2008.
  * [23] Ji Lin, Chuang Gan, and Song Han.  Tsm: Temporal shift module for efficient video understanding.  In Proceedings of ICCV, pages 7083–7093, 2019.
  * [24] Ishan Misra, C Lawrence Zitnick, and Martial Hebert.  Shuffle and learn: unsupervised learning using temporal order verification.  In ECCV, pages 527–544. Springer, 2016.
  * [25] Joe Yue-Hei Ng, Jonghyun Choi, Jan Neumann, and Larry S Davis.  Actionflownet: Learning motion representation for action recognition.  In 2018 WACV, pages 1616–1624. IEEE, 2018.
  * [26] Rui Qian, Tianjian Meng, Boqing Gong, Ming-Hsuan Yang, Huisheng Wang, Serge Belongie, and Yin Cui.  Spatiotemporal contrastive video representation learning.  arXiv preprint arXiv:2008.03800, 2020.
  * [27] Mark Sandler, Andrew Howard, Menglong Zhu, Andrey Zhmoginov, and Liang-Chieh Chen.  Mobilenetv2: Inverted residuals and linear bottlenecks.  In Proceedings of CVPR, pages 4510–4520, 2018.
  * [28] Zheng Shou, Xudong Lin, Yannis Kalantidis, Laura Sevilla-Lara, Marcus Rohrbach, Shih-Fu Chang, and Zhicheng Yan.  Dmc-net: Generating discriminative motion cues for fast compressed video action recognition.  In Proceedings of CVPR, pages 1268–1277, 2019.
  * [29] Karen Simonyan and Andrew Zisserman.  Two-stream convolutional networks for action recognition in videos.  arXiv preprint arXiv:1406.2199, 2014.
  * [30] Du Tran, Lubomir Bourdev, Rob Fergus, Lorenzo Torresani, and Manohar Paluri.  Learning spatiotemporal features with 3d convolutional networks.  In Proceedings of ICCV, pages 4489–4497, 2015.
  * [31] Du Tran, Heng Wang, Lorenzo Torresani, Jamie Ray, Yann LeCun, and Manohar Paluri.  A closer look at spatiotemporal convolutions for action recognition.  In Proceedings of CVPR, pages 6450–6459, 2018.
  * [32] Heng Wang and Cordelia Schmid.  Action recognition with improved trajectories.  In Proceedings of ICCV, pages 3551–3558, 2013.
  * [33] Jiangliu Wang, Jianbo Jiao, Linchao Bao, Shengfeng He, Yunhui Liu, and Wei Liu.  Self-supervised spatio-temporal representation learning for videos by predicting motion and appearance statistics.  In Proceedings of CVPR, pages 4006–4015, 2019.
  * [34] Limin Wang, Yuanjun Xiong, Zhe Wang, Yu Qiao, Dahua Lin, Xiaoou Tang, and Luc Van Gool.  Temporal segment networks: Towards good practices for deep action recognition.  In ECCV, pages 20–36. Springer, 2016.
  * [35] Xiaolong Wang and Abhinav Gupta.  Unsupervised learning of visual representations using videos.  In Proceedings of ICCV, pages 2794–2802, 2015.
  * [36] Saining Xie, Chen Sun, Jonathan Huang, Zhuowen Tu, and Kevin Murphy.  Rethinking spatiotemporal feature learning: Speed-accuracy trade-offs in video classification.  In ECCV, pages 305–321, 2018.
  * [37] Dejing Xu, Jun Xiao, Zhou Zhao, Jian Shao, Di Xie, and Yueting Zhuang.  Self-supervised spatiotemporal learning via video clip order prediction.  In Proceedings of CVPR, pages 10334–10343, 2019.
  * [38] Dejun Zhang, Linchao He, Zhigang Tu, Shifu Zhang, Fei Han, and Boxiong Yang.  Learning motion representation for real-time spatio-temporal action localization.  Pattern Recognition, 103:107312, 2020.

##  Appendix A Appendix

###  A.1 Action to Motion Type Mapping

We have manually annotated all the action classes present in HMDB51 with
motion type classes based on the mapping shown in Table 4. to obtain mHMDB51
dataset. We have used mHMDB51 dataset for training and evaluation of motion
type classifier.

Table 4: Motion Type mapping based on action class Action | Motion Type
---|---
brush_hair | Linear
cartwheel | Projectile
catch | Projectile
chew | Local
clap | Oscillatory
climb | Linear
climb_stairs | Linear
dive | Projectile
draw_sword | Random
dribble | Oscillatory
drink | Local
eat | Local
fall_floor | Random
fencing | Random
flic_flac | Projectile
golf | Projectile
handstand | Projectile
hit | Projectile
hug | Random
jump | Projectile
kick | Random
kick_ball | Random
kiss | Local
laugh | Local
pick | Random
pour | Local
pullup | Oscillatory
punch | Linear
push | Linear
pushup | Oscillatory
ride_bike | Linear
ride_horse | Linear
run | Linear
shake_hands | Local
shoot_ball | Projectile
shoot_bow | Linear
shoot_gun | Local
sit | Random
situp | Oscillatory
smile | Local
smoke | Local
somersault | Projectile
stand | Random
swing_baseball | Projectile
sword | Random
sword_exercise | Random
talk | Local
throw | Projectile
turn | Random
walk | Linear
wave | Local


###### keywords:

Finite elements, Mixed finite elements, MFEM library, Solution comparison,
Laplace problem, Shape functions order, Mesh refinement level

###### keywords:

Elementos finitos, Elementos finitos mixtos, Librería MFEM, Comparación de
soluciones, Problema de Laplace, Refinamiento de malla

[firstpage = 1, volume = 0, number = 0, month = 00, year = 1900, day = 00,
monthreceived = 0, yearreceived = 1900, monthaccepted = 0, yearaccepted =
1900]

] authors

department = Departamento de Matemáticas, institution = Universidad Nacional
de Colombia, city = Bogotá D.C., country = Colombia ]

In this paper, we develop two finite element formulations for the Laplace
problem and find the way in which they are equivalent. Then we compare the
solutions obtained by both formulations, by changing the order of the shape
functions and the refinement level of the mesh (star with rhomboidal
elements). And, we will give an overview of MFEM library from the LLNL
(Lawrence Livermore National Laboratory), as it is the library used to obtain
the solutions.

En este artículo, desarrollamos dos formulaciones de elementos finitos, la de
Lagrange y la mixta, y encontramos la manera en que son equivalentes. Luego,
comparamos las soluciones obtenidas mediante ambas formulaciones al cambiar el
grado de las "shape functions" y el nivel de refinamiento de la malla (una
estrella con elementos romboidales). Y, daremos una revisión general de la
librería MFEM, ya que es la librería utilizada para obtener las soluciones.

65N30

Note: This work was done during the second period of 2020 in the course
"Beyond Research" from the National University of Colombia. It was supervised
by Juan Galvis and Boyan Lazarov.

##  1 Theoretical framework

In this section we are going to study the theoretic background of the project.
First, we are going to review the two finite element methods used (with the
problem they solve) and then, give some information about the library. In the
finite element parts we’ll develop a problem and define the finite element
spaces used; all this in two dimensions. And, for the library part, we’ll give
an overview of its characteristics and the general structure of the code.

###  1.1 Lagrange finite elements

For this method, we consider the following problem [1]:

$\begin{split}-\Delta&p=f\text{ in }\Omega\\\ &p=0\text{ in
}\Gamma\end{split}$ (1)

where $\Omega\subseteq\mathbb{R}^{2}$ is an open-bounded domain with boundary
$\Gamma$, $f$ is a given function and $\Delta p=\frac{\partial^{2}p}{\partial
x^{2}}+\frac{\partial^{2}p}{\partial y^{2}}$. Consider the space $V$:

$V=\\{v:v\text{ continuous on }\Omega,\frac{\partial v}{\partial
x},\frac{\partial v}{\partial y}\text{ piecewise continuous on }\Omega\text{
and }v=0\text{ on }\Gamma\\}$

Now, we can multiply in the first equation of (1) by some $v\in V$ ($v$ is
called test function) and integrate over $\Omega$:

$-\int_{\Omega}\Delta p\ v=\int_{\Omega}f\ v$ (2)

Applying divergence theorem, the following Green’s formula can be deduced [1]:

$-\int_{\Omega}\Delta p\ v=\int_{\Omega}\nabla v\cdot\nabla p-\int_{\Gamma}v\
\nabla p\cdot\eta$ (3)

where $\eta$ is the outward unit normal to $\Gamma$.

Since $v=0$ on $\Gamma$, the third integral equals $0$.

Remark: The boundary integral does not depend on $p$’s value on $\Gamma$ but
rather on it’s derivative in $\Gamma$. And, this is what’s called an essential
boundary condition.

Then, replacing (3) on (2), we get:

$\int_{\Omega}\nabla v\cdot\nabla p=\int_{\Omega}f\ v$ (4)

Note:[1] If $p\in V$ satisfies (4) for all $v\in V$ and is sufficiently
regular, then $p$ also satisfies (1), ie, it’s a solution for our problem.

In order to set the problem for a computer to solve it, we are going to
discretize it and encode it into a linear system.

First, consider a triangulation $T_{h}$ of the domain $\Omega$. This is,
$T_{h}=\\{K_{1},\dots,K_{m}\\}$ a set of non-overlapping triangles such that
$\Omega=K_{1}\cup\dots\cup K_{m}$ and no vertex ($N_{i}$) of one triangle lies
on the edge of another triangle:

Figure 1: Triangulation of $\Omega$

Note: Triangles have been separated in the edges to take a better look, but
the triangulation has no empty spaces.

The $h$ in the notation $T_{h}$ is important for the project because it gives
a sense of the size for the mesh. It is defined as follows:
$h=\max\\{diam(K):K\in T_{h}\\}$ where $diam(K)=\text{longest side of }K$.

Now, let $V_{h}=\\{v:v\text{ continuous on }\Omega,v|_{K}\text{ linear for
}K\in T_{h},\ v=0\text{ on }\Gamma\\}$.
If we consider the nodes ($N_{1},\dots,N_{M}$) of the triangulation that are
not on the boundary, since $v=0$ there, and define the functions
$\varphi_{j}(N_{i})=\left\\{\begin{array}[]{lcc}1&,\ i=j\\\ \\\ 0&,\
i\not=j\\\ \end{array}\right.$ for $i,j=1,\dots,M$ in a way that
$\varphi_{j}\in V_{h}$:

Figure 2: Function $\varphi_{j}$

With this, $V_{h}=gen\\{\varphi_{i}:i=1,\dots,M\\}$ because, for $v(x)\in
V_{h}$,
$v(x)=\sum_{j=1}^{M}\xi_{j}\varphi_{j}(x),$ with $\xi_{j}=v(N_{j})\ and\
x\in\Omega\cup\Gamma$. So, $V_{h}$ is a finite-dimensional subspace of $V$.
[1]

Then, if $p_{h}\in V_{h}$ satisfies (4) for all $v\in V_{h}$ then, in
particular:

$\int_{\Omega}\nabla p_{h}\cdot\nabla\varphi_{j}=\int_{\Omega}f\ \varphi_{j},\
\ j=1,\dots,M$ (5)

As, $\nabla p_{h}=\sum_{i=1}^{M}\xi_{i}\nabla\varphi_{i}$ with
$\xi_{i}=p_{h}(N_{i})$, replacing on (5) we get:

$\sum_{i=1}^{M}\xi_{i}\int_{\Omega}\nabla\varphi_{i}\cdot\nabla\varphi_{j}=\int_{\Omega}f\
\varphi_{j},\ \ j=1,\dots,M$ (6)

Finally, (6) is a linear system of $M$ equations and $M$ unknowns
($\xi_{1},\dots,\xi_{M}$), which can be written as:

$A\xi=b$ (7)

where $A[i,j]=\int_{\Omega}\nabla\varphi_{i}\cdot\nabla\varphi_{j}$,
$\xi[i]=p_{h}(N_{i})$ and $b[i]=\int_{\Omega}f\ \varphi_{i}$.

In [1], it is shown that (7) has an unique solution and that matrix $A$ has
useful properties for computing with it. Also, we can solve (7) with MFEM
library.

###  1.2 Mixed finite elements

First, let’s define some important spaces, where $\Omega$ is a bounded domain
in $\mathbb{R}^{2}$ and $\Gamma$ its boundary [2]:

$L^{2}(\Omega)=\\{v:\Omega\rightarrow\mathbb{R}\
\Big{|}\int_{\Omega}v^{2}<\infty\\}$ $H^{1}(\Omega)=\\{v\in L^{2}(\Omega)\
\Big{|}\ \frac{\partial v}{\partial x},\frac{\partial v}{\partial y}\in
L^{2}(\Omega)\\}$ $H_{0}^{1}(\Omega)=\\{v\in H^{1}(\Omega)\ |\ v=0\ on\
\Gamma\\}$ $H(div;\Omega)=\\{\mathbf{v}\in L^{2}(\Omega)\times L^{2}(\Omega)\
|\ div(\mathbf{v})\in L^{2}(\Omega)\\}$

As above, let $\Omega\in\mathbb{R}^{2}$ be a bounded domain with boundary
$\Gamma$ and consider the following problem [2]:

$\begin{split}-\Delta&p=f\text{ in }\Omega\\\ &p=0\text{ in
}\Gamma\end{split}$ (1)

where $f\in L^{2}(\Omega)$ and $\Delta p=\frac{\partial^{2}p}{\partial
x^{2}}+\frac{\partial^{2}p}{\partial y^{2}}$.

This problem is the same problem considered in 2.1, but with a special
condition for $f$, and can be reduced to:

$\int_{\Omega}\nabla v\cdot\nabla p=\int_{\Omega}f\ v,\text{ for all $v\in
V$}$

where Dirichlet boundary condition ($p=0\ in\ \Gamma$) is essential.

Remark: The space $V$ can be replaced with $H_{0}^{1}(\Omega)$ as seen in [2].

However, for mixed formulation, boundary won’t be essential but natural:

Let $u=\nabla p$ in $\Omega$.
With this, problem (1) can be written as:

$\begin{split}&u=\nabla p\text{ in }\Omega\\\ &div(u)=-f\text{ in }\Omega\\\
&p=0\text{ in }\Gamma\end{split}$ (2)

because $\Delta p=div(\nabla p)$.

Now, following a similar procedure as in section 2.1:

Multiply the first equation of (2) by some $\mathbf{v}\in H(div;\Omega)$ and
integrate both sides:

$\int_{\Omega}u\ \mathbf{v}=\int_{\Omega}\nabla p\cdot\mathbf{v}$ (3)

Consider Green’s identity [2]:

$\int_{\Omega}\mathbf{v}\cdot\nabla p+\int_{\Omega}p\
div(\mathbf{v})=\int_{\Gamma}(\mathbf{v}\cdot\eta)p$ (4)

Replacing (4) in (3), and considering the third equation of (2), we get:

$\int_{\Omega}u\ \mathbf{v}+\int_{\Omega}p\
div(\mathbf{v})=\int_{\Gamma}(\mathbf{v}\cdot\eta)p$ (5)

where $\eta$ is the normal vector exterior to $\Gamma$.

On the other hand, we can multiply the second equation of (2) by some $w\in
L^{2}(\Omega)$, integrate and obtain:

$\int_{\Omega}w\ div(u)=-\int_{\Omega}f\ w$ (6)

Remark: The boundary integral depends directly on the value of $p$ in
$\Gamma$. And, this is what’s called a natural boundary condition.

Finally, applying boundary condition $p=0\ \text{ in }\Gamma$ into (5), and
joining (5) and (6). We get the following problem deduced from (2):

$\begin{split}&\int_{\Omega}u\ \mathbf{v}+\int_{\Omega}p\ div(\mathbf{v})=0\\\
&\int_{\Omega}w\ div(u)=-\int_{\Omega}f\ w\end{split}$ (7)

Note: For this problem, the objective is to find $(u,p)\in H(div;\Omega)\times
L^{2}(\Omega)$ such that (7) is satisfied for all $\mathbf{v}\in
H(div;\Omega),w\in L^{2}(\Omega)$.

For the discretized problem related to (7), define [2] the following spaces
for a fixed triangulation $T_{h}$ of the domain $\Omega$ and a fixed integer
$k\geq 0$:

$\begin{split}&H_{h}:=\\{\mathbf{v_{h}}\in
H(div;\Omega):\mathbf{v_{h}}|_{K}\in RT_{k}(K)\text{ for all }K\in T_{h}\\}\\\
&L_{h}:=\\{w_{h}\in L^{2}(\Omega):w_{h}|_{K}\in\mathbb{P}_{k}(K)\text{ for all
}K\in T_{h}\\}\end{split}$

where

$\begin{split}&\mathbb{P}_{k}(K)=\\{p:K\rightarrow\mathbb{R}\ :\ p\text{ is a
polynomial of degree }\leq k\\}\\\
&RT_{k}(K)=[\mathbb{P}_{k}(K)\times\mathbb{P}_{k}(K)]+\mathbb{P}_{k}(K)x\end{split}$

Note that $\mathbf{p}\in RT_{k}(K)$ if and only if there exist
$p_{0},p_{1},p_{2}\in\mathbb{P}_{k}(K)$ such that

$\mathbf{p}(x)=\begin{pmatrix}p_{1}(x)\\\
p_{2}(x)\end{pmatrix}+p_{0}(x)\begin{pmatrix}x\\\ y\end{pmatrix}\text{ for all
}\begin{pmatrix}x\\\ y\end{pmatrix}\in K$

Also, $\mathbf{p}$ has a degree of $k+1$.

Then, problem (7) can be changed to: find $(u_{h},p_{h})\in H_{h}\times L_{h}$
such that

$\begin{split}&\int_{\Omega}u_{h}\ \mathbf{v}_{h}+\int_{\Omega}p_{h}\
div(\mathbf{v}_{h})=0\\\ &\int_{\Omega}w_{h}\ div(u_{h})=-\int_{\Omega}f\
w_{h}\end{split}$ (8)

for all $\mathbf{v}_{h}\in H_{h},w_{h}\in L_{h}$.

As spaces $H_{h}$ and $L_{h}$ are finite dimensional, they have a finite
basis. That is, $H_{h}=gen\\{\varphi_{i}:i=1,\dots,M\\}$ and
$L_{h}=gen\\{\psi_{j}:j=1,\dots,N\\}$. Then,
$u_{h}=\sum_{i=i}^{M}u_{i}\varphi_{i}$ and
$p_{h}=\sum_{j=1}^{N}p_{j}\psi_{j}$, where $u_{i}$ and $p_{j}$ are scalars.

In particular, as $\varphi_{k}\in H_{h}$ and $\psi_{l}\in L_{h}$, we have that
problem (8) can be written as

$\begin{split}&\int_{\Omega}\left(\sum_{i=i}^{M}u_{i}\varphi_{i}\right)\varphi_{k}+\int_{\Omega}\left(\sum_{j=1}^{N}p_{j}\psi_{j}\right)div(\varphi_{k})=0\\\
&\int_{\Omega}\psi_{l}div\left(\sum_{i=1}^{M}u_{i}\varphi_{i}\right)=\int_{\Omega}f\psi_{l}\
\end{split}$ (9)

for $k=1,\dots,M$ and $l=1,\dots,N$. Which is equivalent to the following by
rearranging scalars:

$\begin{split}&\sum_{i=i}^{M}u_{i}\int_{\Omega}\varphi_{i}\cdot\varphi_{k}+\sum_{j=1}^{N}p_{j}\int_{\Omega}\psi_{j}div(\varphi_{k})=0\\\
&\sum_{i=i}^{M}u_{i}\int_{\Omega}\psi_{l}div(\varphi_{i})=\int_{\Omega}f\psi_{l}\end{split}$
(10)

for $k=1,\dots,M$ and $l=1,\dots,N$. This problem (10) can be formulated into
the following matrix system

$\begin{pmatrix}A&B\\\ B^{t}&0\end{pmatrix}\begin{pmatrix}U\\\
P\end{pmatrix}=\begin{pmatrix}0\\\ F\end{pmatrix}$ (11)

where $A$ is a $N\times N$ matrix, $B$ is a $M\times N$ matrix with $B^{t}$
it’s transpose, $U$ is a $M$-dimensional column vector and $P,F$ are
$N$-dimensional column vectors.
The entries of these arrays are
$A[i,j]=\int_{\Omega}\varphi_{i}\cdot\varphi_{j}$,
$B[i,j]=\int_{\Omega}\psi_{j}div(\varphi_{i})$, $U[i]=u_{i}$, $P[i]=p_{i}$ and
$F[i]=\int_{\Omega}f\psi_{i}$.

(11) is a multilinear system that can be solved for $(U,P)$ with a computer
using MFEM library. Note that with the entries of $U$ and $P$, the solution
$(u_{h},p_{h})$ of (8) can be computed by their basis representation.

Note: The spaces defined to discretize the problem are called Raviart-Thomas
finite element spaces. The fixed integer k is also called the order of the
shape functions. And, the parameter $h$ is the same as in section 2.1, which
is a meassure of size for $T_{h}$.

###  1.3 Finite elements summary

In sections 2.1 and 2.2 we studied two finite element methods. In general
aspects, this is what was done:

  * •

Consider the problem of solving Poisson’s equation with homogeneous Dirichlet
boundary conditions. That is, the problem considered in previous sections.

  * •

Multiply by some function (test function) and integrate.

  * •

Develop some equations applying boundary conditions.

  * •

Discretize the domain.

  * •

Define some finite-dimensional function spaces.

  * •

Assemble the basis into the equation and form a matrix system.

The functions that form part of the finite-dimensional spaces are called
$shape\ functions$. In Lagrange formulation, those where the functions in
$V_{h}$, and in mixed formulation, those where the functions in $H_{h}$ and
$L_{h}$.

The parameter $h$, denotes the size of the elements in the triangulation of
the domain.

Both problems were solved with Dirichlet boundary condition ($=0$). In
Lagrange formulation it was essential, and in mixed formulation, it was
natural.

In a more general aspect, the discretization of the space can be done without
using triangles, but rather using quads or other figures.

###  1.4 Higher order shape functions

This is a very brief section that has the purpose of explaining a little bit
of finite elements order, because in section 3 we will use different orders
for the shape functions.

In general aspects, the order of a shape function is similar to the order of a
polynomial. In mixed formulation we approached this when talking about
Raviart-Thomas spaces, as in this spaces if the order of the polynomial is
$k$, then the order of the shape function is $k+1$.

In the original introduction of the Lagrange formulation, the order of the
shape functions was set to one. Better approximations can be obtained by using
polynomials of higher order. Instead of defining

$V_{h}=\\{v:v\text{ continuous on }\Omega,v|_{K}\text{ linear for }K\in
T_{h},\ v=0\text{ on }\Gamma\\}$

one can define, for a fixed order $k$:

$V^{k}_{h}=\\{v:v\text{ continuous on }\Omega,v|_{K}\text{ polynomial of order
at most }K\in T_{h},\ v=0\text{ on }\Gamma\\}.$

Remark: For a fixed $k$, Lagrange shape functions have order 1 less than mixed
shape functions.

For example, as seen in [3], the space of Bell triangular finite elements for
a given triangulation $T_{h}$ is the space of functions that are polynomials
of order 5 when restricted to every triangle $K\in T_{h}$. That is, if $v$ is
in this space, then:

$v|_{K}(x,y)=a_{1}x^{5}+a_{2}y^{5}+a_{3}x^{4}y+a_{4}xy^{4}+\dots+a_{16}x+a_{17}y+a_{18}$

for all $K\in T_{h}$. Here, the constants $a_{i},\ i=1,\dots,18$ correspond to
$v$’s DOF (degrees of freedom).

Figure 3: Finite element of order 2

Figure 4: Finite elements of orders 5 (left) and 10 (right)

###  1.5 MFEM library

In this project, we worked with MFEM’s Example#1 and Example#5 which can be
found on [4]. Example#1 uses standard Lagrange finite elements and Example#5
uses Raviart-Thomas mixed finite elements. Further, in section 3.1, we find
the parameters so that both problems are equivalent and then (section 3.4), we
compare the solutions.

####  1.5.1 Overview

According to it’s official site [4], MFEM is a free, lightweight, scalable C++
library for finite element methods that can work with arbitrary high-order
finite element meshes and spaces.

MFEM has a serial version (which we are using) and a parallel version (for
parallel computation).

The main classes (with a brief and superficial explanation of them) that we
are going to use in the code are:

  * •

Mesh: domain with the partition.

  * •

FiniteElementSpace: space of functions defined on the finite element mesh.

  * •

GridFunction: mesh with values (solutions).

  * •

$\\_$Coefficient: values of GridFunctions or constants.

  * •

LinearForm: maps an input function to a vector for the rhs.

  * •

BilinearForm: used to create a global sparse finite element matrix for the
lhs.

  * •

$\\_$Vector: vector.

  * •

$\\_$Solver: algorithm for solution calculation.

  * •

$\\_$Integrator: evaluates the bilinear form on element’s level.

The ones that have $\\_$ are various classes whose name ends up the same and
work similarly.

Note:
lhs: left hand side of the linear system.
rhs: right hand side of the linear system.

####  1.5.2 Code structure

An MFEM general code has the following steps (directly related classes with
the step are written):

  1. 1.

Receive archive (.msh) input with the mesh and establish the order for the
finite element spaces.

  2. 2.

Create mesh object, get the dimension, and refine the mesh (refinement is
optional). Mesh

  3. 3.

Define the finite element spaces required. FiniteElementSpace

  4. 4.

Define coefficients, functions, and boundary conditions of the problem.
XCoefficient

  5. 5.

Define the LinearForm for the rhs and assemble it. LinearForm, XIntegrator

  6. 6.

Define the BilinearForm for the lhs and assemble it. BilinearForm, XIntegrator

  7. 7.

Solve the linear system. XSolver, XVector

  8. 8.

Recover solution. GridFunction

  9. 9.

Show solution with a finite element visualization tool like Glvis (optional).

##  2 A case study

In this section: we take examples 1 and 5 from [4], define their problem
parameters in such way that they’re equivalent, create a code that implements
both of them at the same time and compares both solutions ($L_{2}$ norm), run
the code with different orders, and analyse the results.

Some considerations to have into account are:

  * •

For a fair comparison, order for Mixed method should be 1 less than order for
Lagrange method. Because, with this, both shape functions would have the same
degree.

  * •

The code has more steps than shown in section 2.3.2 because we are running two
methods and comparing solutions.

  * •

We will compare pressures and velocities with respect to the order of the
shape functions and the size of the mesh ($h$ parameter).

  * •

For the problem, the exact solution is known, so, we will use it for
comparison.

  * •

The max order and refinement level to be tested is determined by our
computational capacity (as long as solvers converge fast).

  * •

The mesh used is a star with rhomboidal elements.

###  2.1 Problem

Example#1 [4]:

$\begin{split}-\Delta&p=1\text{ in }\Omega\\\ &p=0\text{ in
}\Gamma\end{split}$ (1)

Example#5 [4]:

$\begin{split}&k\mathbf{u}+\nabla p=f\text{ in }\Omega\\\
&-div(\mathbf{u})=g\text{ in }\Omega\\\ &-p=p_{0}\text{ in }\Gamma\end{split}$
(2)

From the first equation of (2):

$\mathbf{u}=\frac{f-\nabla p}{k}$ (3)

Then, replacing (3) on the second equation of (2):

$-div\left(\frac{f-\nabla p}{k}\right)=g$ (4)

If we set $k=1;\ f=0\ and\ g=-1$ in (4), we get:

$-\Delta p=1$ (5)

which is the first equation of (1).

So, setting ($*$) $p_{0}=0,\ k=1;\ f=0\ and\ g=-1$ in (2), we get:

$\begin{split}&\mathbf{u}+\nabla p=0\text{ in }\Omega\\\
&-div(\mathbf{u})=-1\text{ in }\Omega\\\ &-p=0\text{ in }\Gamma\end{split}$
(6)

Notice that from the first equation we get that $\mathbf{u}=-\nabla p$. This
is important because in problem (1) we don’t get $\mathbf{u}$ solution from
the method, so, in the code, we will have to find it from $p$’s derivatives.

In the code, we will set the value of the parameters in the way shown here, so
that both problems are the same. As seen in (3)-(5), problem (6) is equivalent
to problem (1) with the values assigned for coefficients and functions in
($*$).

###  2.2 Code

The first part of the code follows the structure mentioned in 2.3.2, but
implemented for two methods at the same time (and with some extra lines for
comparison purposes). Also, when defining boundary conditions, the essential
one is established different from the natural one. And, after getting all the
solutions, there’s a second part of the code where solutions are compared
between them and with the exact one.

Note:
The complete code with explanations can be found on the Appendix A.

However, before taking a look into it, here’s the convention used for
important variable names along the code:

Notation:

Variable Name | Object
---|---
X_space | Finite element space X
X_mixed | Variable assigned to a mixed method related object
u | Velocity solution
p | Pressure solution
X_ex | Variable assigned to an exact solution object

###  2.3 Tests

The tests will be run on the following domain:

Figure 5: Star domain for tests

Each run test is determined by the order of Lagrange shape functions and the h
parameter of the mesh. Remember that mixed shape functions have order equal to
$\textit{order}-1$. The parameter order is changed directly from the command
line, while the parameter h is changed via the number of times that the mesh
is refined ($h=h(\\#refinements)$). As we refine the mesh more times, finite
elements of the partition decrease their size, and so, the parameter $h$
decreases.

Tests will be made with: $order=1,\dots,N$ and $refinements=0,\dots,M$, where
$N,M$ depend on the computation capacity. The star mesh comes with a default
partition which is shown below:

Figure 6: Mesh with no refinement

Results will be presented in graphs. However, all the exact values that were
computed can be found in the Appendix B.

###  2.4 Results

Before showing the graphs, this is the output received in the visualization
tool (Glvis) when running the code with $\textit{order}=2$ and
$\\#Refinements=3$ (graphically, Lagrange and Mixed solutions look the same):

Figure 7: Glvis Visualization: Pressure (left) and Velocity (right)

Note: Although velocity is a vector on each point, Glvis visualization tool
doesn’t shows it like that. It rather shows the $L^{2}$ norm of the vector.

In the following graphs, if $u=(u_{x},u_{y})$ is the solution obtained by the
mixed or Lagrange finite element method and $u_{ex}=(u_{x_{ex}},u_{y_{ex}})$
is the exact solution for the problem, then:

$U_{error}=\frac{\sqrt{\left(||u_{x}-u_{x_{ex}}||_{L^{2}}\right)^{2}+\left(||u_{y}-u_{y_{ex}}||_{L^{2}}\right)^{2}}}{||u_{ex}||_{L^{2}}}$

Figure 8: Order = 1 Figure 9: Order = 2

Figure 10: Order = 3 Figure 11: Order = 4

###  2.5 Analysis

This section was done by analyzing the tables presented on the Appendix B.

To understand the information presented, take into account the following:

  * •

The exact solution would have value $1$ in X err.

  * •

If the two solutions obtained (Lagrange and Mixed) are exactly the same, the
value in P comp and U comp would be $0$.

  * •

Lower values of $h$ mean more mesh refinements, ie, smaller partition
elements.

As it was expected, computational time increases as order and refinements
increase.

Here are the most relevant observations that can be obtained after analysing
the data corresponding to absolute errors:

  * •

For fixed order, absolute errors have little variation when reducing $h$ (max
variation is $4.722$e$-03$ in $Uerr$ order 1).

  * •

Absolute errors variation (respect to refinement) is lower when order is
higher. For example; in order 2, $Perr$ is the same for each $h$ (up tu three
decimal places); while in order 6, $Perr$ is the same for each $h$ (up to five
decimal places).

  * •

For fixed $h$, absolute errors remain almost constant between orders.

  * •

$Perr$ (absolute error obtained for pressure with Lagrange) is always lower
than $Pmx\ err$ (absolute error obtained for pressure with mixed).

  * •

For fixed order, $Perr$ increases as $h$ decreases, while $Pmx\ err$ decreases
as $h$ decreases.

  * •

$Uerr$ (absolute error obtained for velocity with Lagrange) is always lower
than $Umx\ err$ (absolute error obtained for velocity with mixed).

  * •

For fixed order, $Uerr$ increases as $h$ decreases, while $Umx\ err$ decreases
as $h$ decreases.

  * •

As order increases, pressure absolute errors tend to be the same. In order 10,
the difference between $Perr$ and $Pmx\ err$ is $0.000001$.

  * •

As order increases, velocity absolute errors tend to be the same. In order 10,
the difference between $Uerr$ and $Umx\ err$ is $<0.0000009$.

And now, the most relevant observations that can be obtained after analysing
the data corresponding to comparison errors:

  * •

Comparison errors, $Ucomp$ and $Pcomp$, decrease as $h$ decreases.

  * •

When order increases, comparisons errors are lower for fixed $h$.

  * •

Comparison error tends to $0$, as expected.

  * •

Pressure comparison error lowers faster than velocity comparison error.
Maximum comparison errors were found on order 1 with no refinements, where
$Pcomp\approx 7.5$e$-02$ and $Ucomp\approx 3.7$e$-02$, and in minimum
comparison errors were found on order 10 with 1 refinement (higher refinement
level computed for order 10), where $Pcomp\approx 5.1$e$-06$ and $Ucomp\approx
9.8$e$-04$. It can be seen that $Pcomp$ improved in almost four decimal places
while $Ucomp$ improved in just 2.

  * •

For a fixed order, comparison error can be similar to a higher order
comparison error, as long as enough refinements are made.

##  3 Conclusion

Adding up to the observations made in section 3.5, Lagrange solution and mixed
solution tend to be the same when order and refinement levels increase, as
expected. Also, Lagrange formulation is implemented more easily compared to
mixed formulation but, with mixed formulation one can obtain pressure and
velocity solutions at once. Furthermore, in MFEM, natural boundary conditions
can be forced in an easier way compared to essential boundary conditions.
Finally, it’s important to note that finite element methods are a powerful
mathematical tool used to solve potentially difficult problems.

## References

  * [1] Claes Johnson. Numerical Solution of Partial Differential Equations by the Finite Element Method. ISBN10 048646900X. Dover Publications Inc. 2009.
  * [2] Gabriel N. Gatica. A Simple Introduction to the Mixed Finite Element Method. Theory and Applications. ISBN 978-3-319-03694-6. Springer. 2014.
  * [3] Juan Galvis & Henrique Versieux. Introdução à Aproximação Numérica de Equações Diferenciais Parciais Via o Método de Elementos Finitos. ISBN: 978-85-244-325-5. 28 Colóquio Brasileiro de Matemática. 2011.
  * [4] MFEM. Principal online page at: mfem.org. Code Documentation. Examples #1 and #5.

##  4 Appendices

###  4.1 Appendix A

Here, the code used (written in C++) is shown, with a brief explanations of
it’s functionality.

$\triangleright$Include the required libraries (including MFEM) and begin main
function.

⬇

#include "mfem.hpp"

#include <fstream>

#include <iostream>

using namespace std;

using namespace mfem;

int main(int argc, char *argv[]){

$\triangleright$Parse command-line options (in this project we only change
"order" option) and print them.

⬇

const char *mesh_file = "../data/star.mesh";

int order = 1;

bool visualization = true;

OptionsParser args(argc, argv);

args.AddOption(&mesh_file, "-m", "–mesh",

"Mesh␣file␣to␣use.");

args.AddOption(&order, "-o", "–order",

"Finite␣element␣order␣(polynomial␣degree).");

args.AddOption(&visualization, "-vis", "–visualization", "-no-vis", "–no-
visualization",

"Enable␣or␣disable␣GLVis␣visualization.");

args.Parse();

if (!args.Good()){

args.PrintUsage(cout);

return 1;

}

args.PrintOptions(cout);

$\triangleright$Create mesh object from the star.mesh archive and get it’s
dimension.

⬇

Mesh *mesh = new Mesh(mesh_file,1,1);

int dim = mesh->Dimension();

$\triangleright$Refine the mesh a given number of times (uniform refinement).

⬇

int ref_levels;

cout << "Refinements:␣";

cin >> ref_levels;

for (int l = 0; l < ref_levels; l++){

mesh->UniformRefinement();

}

$\triangleright$Get size indicator for mesh size (h_max) and print it.

⬇

double mesh_size, h = 0;

for (int i=0;i<mesh->GetNE();i++){

mesh_size = mesh->GetElementSize(i,2);

if(mesh_size>h){

h = mesh_size;

}

}

cout << "h:␣" << h << endl;

$\triangleright$Define finite element spaces. For mixed finite element method,
the order will be one less than for Lagrange finite element method. The last
one is a vector L2 space that we will use later to get mixed velocity
components.

⬇

FiniteElementCollection *H1 = new H1_FECollection(order,dim);

FiniteElementSpace *H1_space = new FiniteElementSpace(mesh,H1);

FiniteElementCollection *hd(new RT_FECollection(order-1,dim));

FiniteElementCollection *l2(new L2_FECollection(order-1,dim));

FiniteElementSpace *Hdiv_space = new FiniteElementSpace(mesh,hd);

FiniteElementSpace *L2_space = new FiniteElementSpace(mesh,l2);

FiniteElementSpace *V_space = new FiniteElementSpace(mesh,l2,2);

$\triangleright$Define the parameters of the mixed problem. C functions are
defined at the end. Boundary condition is natural.

⬇

ConstantCoefficient k(1.0);

void fFun(const Vector & x, Vector & f);

VectorFunctionCoefficient fcoeff(dim, fFun);

double gFun(const Vector & x);

FunctionCoefficient gcoeff(gFun);

double f_bound(const Vector & x);

FunctionCoefficient fbndcoeff(f_bound);

$\triangleright$Define the parameters of the Lagrange problem. Boundary
condition is essential.

⬇

ConstantCoefficient one(1.0);

Array<int> ess_tdof_list;

if (mesh->bdr_attributes.Size()){

Array<int> ess_bdr(mesh->bdr_attributes.Max());

ess_bdr = 1;

H1_space->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);

}

$\triangleright$Define the exact solution. C functions are defined at the end.

⬇

void u_ex(const Vector & x, Vector & u);

double p_ex(const Vector & x);

double u_ex_x(const Vector & x);

double u_ex_y(const Vector & x);

$\triangleright$Get space dimensions and crate vectors for the right hand
side.

⬇

Array<int> block_offsets(3);

block_offsets[0] = 0;

block_offsets[1] = Hdiv_space->GetVSize();

block_offsets[2] = L2_space->GetVSize();

block_offsets.PartialSum();

BlockVector rhs_mixed(block_offsets);

Vector rhs(H1_space->GetVSize());

$\triangleright$Define the right hand side. These are LinearForm objects
associated to some finite element space and rhs vector. "f" and "g" are for
the mixed method and "b" is for the other method. "rhs" vectors are the
variables that store the information of the right hand side.

⬇

LinearForm *fform(new LinearForm);

fform->Update(Hdiv_space, rhs_mixed.GetBlock(0), 0);

fform->AddDomainIntegrator(new VectorFEDomainLFIntegrator(fcoeff));

fform->AddBoundaryIntegrator(new VectorFEBoundaryFluxLFIntegrator(fbndcoeff));

fform->Assemble();

LinearForm *gform(new LinearForm);

gform->Update(L2_space, rhs_mixed.GetBlock(1), 0);

gform->AddDomainIntegrator(new DomainLFIntegrator(gcoeff));

gform->Assemble();

LinearForm *b(new LinearForm);

b->Update(H1_space, rhs, 0);

b->AddDomainIntegrator(new DomainLFIntegrator(one));

b->Assemble();

$\triangleright$Create variables to store the solution. "x" is the vector used
as input in the iterative method.

⬇

BlockVector x_mixed(block_offsets);

GridFunction u_mixed(Hdiv_space), p_mixed(L2_space), ux_mixed(L2_space),
uy_mixed(L2_space), ue(V_space);

Vector x(H1_space->GetVSize());

GridFunction ux(L2_space),uy(L2_space),p(H1_space);

$\triangleright$Define the left hand side for mixed method. This is the
bilinear form representing the Darcy matrix. VectorFEMMassIntegrator is
asociated to $k*u-\nabla p$ and VectorFEDDivergenceIntegrator is asociated to
$div(u)$.

⬇

BilinearForm *mVarf(new BilinearForm(Hdiv_space));

MixedBilinearForm *bVarf(new MixedBilinearForm(Hdiv_space, L2_space));

mVarf->AddDomainIntegrator(new VectorFEMassIntegrator(k));

mVarf->Assemble();

mVarf->Finalize();

SparseMatrix &M(mVarf->SpMat());

bVarf->AddDomainIntegrator(new VectorFEDivergenceIntegrator);

bVarf->Assemble();

bVarf->Finalize();

SparseMatrix & B(bVarf->SpMat());

B *= -1.;

SparseMatrix *BT = Transpose(B);

BlockMatrix D(block_offsets);

D.SetBlock(0,0, &M);

D.SetBlock(0,1, BT);

D.SetBlock(1,0, &B);

$\triangleright$Define the left hand side for Lagrange method. This is the
bilinear form asociated to the laplacian operator. DiffusionIntegrator is
asociated to $\Delta u$. The method FormLinearSystem is only used to establish
the essential boundary condition.

⬇

OperatorPtr A;

Vector XX,BB;

BilinearForm *a(new BilinearForm(H1_space));

a->AddDomainIntegrator(new DiffusionIntegrator(one));

a->Assemble();

a->FormLinearSystem(ess_tdof_list, p, *b, A, XX, BB);

$\triangleright$Solve linear systems with MINRES (for mixed) and CG (for
Lagrange). SetOperator method establishes the lhs. Mult method executes the
iterative algorithm and receives as input: the rhs and the vector to store the
solution. Then convergence result is printed.

⬇

int maxIter(10000);

double rtol(1.e-6);

double atol(1.e-10);

MINRESSolver Msolver;

Msolver.SetAbsTol(atol);

Msolver.SetRelTol(rtol);

Msolver.SetMaxIter(maxIter);

Msolver.SetPrintLevel(0);

Msolver.SetOperator(D);

x_mixed = 0.0;

Msolver.Mult(rhs_mixed, x_mixed);

if (Msolver.GetConverged())

std::cout << "MINRES␣converged␣in␣" << Msolver.GetNumIterations() <<
"␣iterations␣with␣a␣residual␣norm␣of␣" << Msolver.GetFinalNorm() << ".\n";

else

std::cout << "MINRES␣did␣not␣converge␣in␣" << Msolver.GetNumIterations() <<
"␣iterations.␣Residual␣norm␣is␣" << Msolver.GetFinalNorm() << ".\n";

CGSolver Lsolver;

Lsolver.SetAbsTol(atol);

Lsolver.SetRelTol(rtol);

Lsolver.SetMaxIter(maxIter);

Lsolver.SetPrintLevel(0);

Lsolver.SetOperator(*A);

x = 0.0;

Lsolver.Mult(rhs,x);

if (Lsolver.GetConverged())

std::cout << "CG␣converged␣in␣" << Lsolver.GetNumIterations() <<
"␣iterations␣with␣a␣residual␣norm␣of␣" << Lsolver.GetFinalNorm() << ".\n";

else

std::cout << "CG␣did␣not␣converge␣in␣" << Lsolver.GetNumIterations() <<
"␣iterations.␣Residual␣norm␣is␣" << Lsolver.GetFinalNorm() << ".\n";

$\triangleright$Save the solution into GridFunctions, which are used for error
computation and visualization.

⬇

u_mixed.MakeRef(Hdiv_space, x_mixed.GetBlock(0), 0);

p_mixed.MakeRef(L2_space, x_mixed.GetBlock(1), 0);

p.MakeRef(H1_space,x,0);

$\triangleright$Get missing velocities from the solutions obtained. Remember
that $u=-\nabla p$. Mixed components are extracted using the auxiliary
variable "ue" defined before.

⬇

p.GetDerivative(1,0,ux);

p.GetDerivative(1,1,uy);

ux *= -1;

uy *= -1;

VectorGridFunctionCoefficient uc(&u_mixed);

ue.ProjectCoefficient(uc);

GridFunctionCoefficient ux_mixed_coeff(&ue,1);

GridFunctionCoefficient uy_mixed_coeff(&ue,2);

ux_mixed.ProjectCoefficient(ux_mixed_coeff);

uy_mixed.ProjectCoefficient(uy_mixed_coeff);

$\triangleright$Create the asociated Coefficient objects for error
computation.

⬇

GridFunction* pp = &p;

GridFunctionCoefficient p_coeff(pp);

GridFunction* uxp = &ux;

GridFunction* uyp = &uy;

GridFunctionCoefficient ux_coeff(uxp);

GridFunctionCoefficient uy_coeff(uyp);

FunctionCoefficient pex_coeff(p_ex);

VectorFunctionCoefficient uex_coeff(dim,u_ex);

FunctionCoefficient uex_x_coeff(u_ex_x);

FunctionCoefficient uex_y_coeff(u_ex_y);

$\triangleright$Define integration rule.

⬇

int order_quad = max(2, 2*order+1);

const IntegrationRule *irs[Geometry::NumGeom];

for (int i=0; i < Geometry::NumGeom; ++i){

irs[i] = &(IntRules.Get(i, order_quad));

}

$\triangleright$Compute exact solution norms.

⬇

double norm_p = ComputeLpNorm(2., pex_coeff, *mesh, irs);

double norm_u = ComputeLpNorm(2., uex_coeff, *mesh, irs);

double norm_ux = ComputeLpNorm(2., uex_x_coeff, *mesh, irs);

double norm_uy = ComputeLpNorm(2., uex_y_coeff, *mesh, irs);

$\triangleright$Compute absolute errors and print them.

⬇

double abs_err_u_mixed = u_mixed.ComputeL2Error(uex_coeff,irs);

printf("Velocity␣Mixed␣Absolute␣Error:␣%e\n", abs_err_u_mixed / norm_u);

double abs_err_p_mixed = p_mixed.ComputeL2Error(pex_coeff,irs);

printf("Pressure␣Mixed␣Absolute␣Error:␣%e\n", abs_err_p_mixed / norm_p);

double abs_err_p = p.ComputeL2Error(pex_coeff,irs);

printf("Pressure␣Absolute␣Error:␣%e\n", abs_err_p / norm_p);

double abs_err_ux = ux.ComputeL2Error(uex_x_coeff,irs);

double abs_err_uy = uy.ComputeL2Error(uex_y_coeff,irs);

double abs_err_u = pow(pow(abs_err_ux,2)+pow(abs_err_uy,2),0.5);

printf("Velocity␣Absolute␣Error:␣%e\n", abs_err_u / norm_u);

$\triangleright$Compute and print comparison errors.

⬇

double err_ux = ux_mixed.ComputeL2Error(ux_coeff,irs);

double err_uy = uy_mixed.ComputeL2Error(uy_coeff,irs);

double err_u = pow(pow(err_ux,2)+pow(err_uy,2),0.5);

printf("Velocity␣Comparison␣Error:␣%e\n", err_u / norm_u);

double err_p = p_mixed.ComputeL2Error(p_coeff, irs);

printf("Pressure␣Comparison␣Error:␣%e\n", err_p / norm_p);

$\triangleright$Visualize the solutions and the domain.

⬇

char vishost[] = "localhost";

int visport = 19916;

if(visualization){

Vector x_domain(H1_space->GetVSize());

GridFunction domain(H1_space);

x_domain=0.0;

domain.MakeRef(H1_space,x_domain,0);

socketstream dom_sock(vishost, visport);

dom_sock.precision(8);

dom_sock << "solution\n" << *mesh << domain << "window_title␣’Domain’" <<
endl;

socketstream um_sock(vishost, visport);

um_sock.precision(8);

um_sock << "solution\n" << *mesh << u_mixed << "window_title␣’Velocity␣Mixed’"
<< endl;

socketstream pm_sock(vishost, visport);

pm_sock.precision(8);

pm_sock << "solution\n" << *mesh << p_mixed << "window_title␣’Pressure␣Mixed’"
<< endl;

socketstream uxm_sock(vishost, visport);

uxm_sock.precision(8);

uxm_sock << "solution\n" << *mesh << ux_mixed <<
"window_title␣’X␣Velocity␣Mixed’" << endl;

socketstream uym_sock(vishost, visport);

uym_sock.precision(8);

uym_sock << "solution\n" << *mesh << uy_mixed <<
"window_title␣’Y␣Velocity␣Mixed’" << endl;

socketstream p_sock(vishost, visport);

p_sock.precision(8);

p_sock << "solution\n" << *mesh << p << "window_title␣’Pressure’" << endl;

socketstream ux_sock(vishost, visport);

ux_sock.precision(8);

ux_sock << "solution\n" << *mesh << ux << "window_title␣’X␣Velocity’" << endl;

socketstream uy_sock(vishost, visport);

uy_sock.precision(8);

uy_sock << "solution\n" << *mesh << uy << "window_title␣’Y␣Velocity’" << endl;

}

}

$\triangleright$Define C functions.

⬇

void fFun(const Vector & x, Vector & f){

f = 0.0;

}

double gFun(const Vector & x){

return -1.0;

}

double f_bound(const Vector & x){

return 0.0;

}

void u_ex(const Vector & x, Vector & u){

double xi(x(0));

double yi(x(1));

double zi(0.0);

u(0) = - exp(xi)*sin(yi)*cos(zi);

u(1) = - exp(xi)*cos(yi)*cos(zi);

}

double u_ex_x(const Vector & x){

double xi(x(0));

double yi(x(1));

double zi(0.0);

return -exp(xi)*sin(yi)*cos(zi);

}

double u_ex_y(const Vector & x){

double xi(x(0));

double yi(x(1));

double zi(0.0);

return -exp(xi)*cos(yi)*cos(zi);

}

double p_ex(const Vector & x){

double xi(x(0));

double yi(x(1));

double zi(0.0);

return exp(xi)*sin(yi)*cos(zi);

}

###  4.2 Appendix B

The order parameter will be fixed for each table and $h$ parameter is shown in
the first column. To interpret the results take into account that P refers to
pressure, U refers to velocity, mx refers to mixed (from mixed finite element
method), err refers to absolute error (compared to the exact solution), and
comp refers to comparison (the error between the two solutions obtained by the
two different methods).

Order = 1

h | P comp | P err | Pmx err | U comp | U err | U mx err
---|---|---|---|---|---|---
0.572063 | 7.549479e-02 | 1.021287e+00 | 1.025477e+00 | 3.680827e-02 | 1.029378e+00 | 1.037635e+00
0.286032 | 3.627089e-02 | 1.022781e+00 | 1.023990e+00 | 1.727281e-02 | 1.032760e+00 | 1.035055e+00
0.143016 | 1.791509e-02 | 1.023236e+00 | 1.023596e+00 | 9.222996e-03 | 1.033725e+00 | 1.034369e+00
0.0715079 | 8.922939e-03 | 1.023372e+00 | 1.023480e+00 | 5.111295e-03 | 1.033999e+00 | 1.034182e+00
0.035754 | 4.455715e-03 | 1.023412e+00 | 1.023445e+00 | 2.859769e-03 | 1.034077e+00 | 1.034130e+00
0.017877 | 2.226845e-03 | 1.023424e+00 | 1.023435e+00 | 1.603788e-03 | 1.034100e+00 | 1.034115e+00

Order = 2

h | P comp | P err | Pmx err | U comp | U err | U mx err
---|---|---|---|---|---|---
0.572063 | 8.069013e-03 | 1.023329e+00 | 1.023554e+00 | 1.399079e-02 | 1.033924e+00 | 1.034255e+00
0.286032 | 2.138257e-03 | 1.023391e+00 | 1.023470e+00 | 7.845012e-03 | 1.034056e+00 | 1.034146e+00
0.143016 | 5.704347e-04 | 1.023417e+00 | 1.023442e+00 | 4.400448e-03 | 1.034093e+00 | 1.034120e+00
0.0715079 | 1.537926e-04 | 1.023426e+00 | 1.023434e+00 | 2.469526e-03 | 1.034104e+00 | 1.034112e+00
0.035754 | 4.194302e-05 | 1.023428e+00 | 1.023431e+00 | 1.385966e-03 | 1.034107e+00 | 1.034110e+00

Order = 3

h | P comp | P err | Pmx err | U comp | U err | U mx err
---|---|---|---|---|---|---
0.572063 | 8.691241e-04 | 1.023389e+00 | 1.023471e+00 | 8.745151e-03 | 1.034060e+00 | 1.034143e+00
0.286032 | 2.477673e-04 | 1.023417e+00 | 1.023443e+00 | 4.911967e-03 | 1.034094e+00 | 1.034120e+00
0.143016 | 7.316263e-05 | 1.023426e+00 | 1.023434e+00 | 2.756849e-03 | 1.034104e+00 | 1.034112e+00
0.0715079 | 2.178864e-05 | 1.023428e+00 | 1.023431e+00 | 1.547232e-03 | 1.034108e+00 | 1.034110e+00

Order = 4

h | P comp | P err | Pmx err | U comp | U err | U mx err
---|---|---|---|---|---|---
0.572063 | 3.199774e-04 | 1.023412e+00 | 1.023448e+00 | 6.119857e-03 | 1.034088e+00 | 1.034124e+00
0.286032 | 9.547574e-05 | 1.023424e+00 | 1.023435e+00 | 3.434952e-03 | 1.034103e+00 | 1.034114e+00
0.143016 | 2.862666e-05 | 1.023428e+00 | 1.023431e+00 | 1.927814e-03 | 1.034107e+00 | 1.034111e+00

Order = 5

h | P comp | P err | Pmx err | U comp | U err | U mx err
---|---|---|---|---|---|---
0.572063 | 1.552006e-04 | 1.023420e+00 | 1.023439e+00 | 4.578518e-03 | 1.034099e+00 | 1.034117e+00
0.286032 | 4.658038e-05 | 1.023427e+00 | 1.023433e+00 | 2.569749e-03 | 1.034106e+00 | 1.034112e+00
0.143016 | 1.406993e-05 | 1.023429e+00 | 1.023431e+00 | 1.442205e-03 | 1.034108e+00 | 1.034110e+00

Order = 6

h | P comp | P err | Pmx err | U comp | U err | U mx err
---|---|---|---|---|---|---
0.572063 | 8.612580e-05 | 1.023424e+00 | 1.023435e+00 | 3.584133e-03 | 1.034103e+00 | 1.034114e+00
0.286032 | 2.600417e-05 | 1.023428e+00 | 1.023431e+00 | 2.011608e-03 | 1.034107e+00 | 1.034111e+00
0.143016 | 7.897631e-06 | 1.023429e+00 | 1.023430e+00 | 1.128989e-03 | 1.034109e+00 | 1.034110e+00

Order = 7

h | P comp | P err | Pmx err | U comp | U err | U mx err
---|---|---|---|---|---|---
0.572063 | 5.243187e-05 | 1.023426e+00 | 1.023433e+00 | 2.899307e-03 | 1.034105e+00 | 1.034112e+00
0.286032 | 1.589631e-05 | 1.023429e+00 | 1.023431e+00 | 1.627221e-03 | 1.034108e+00 | 1.034110e+00

Order = 8

h | P comp | P err | Pmx err | U comp | U err | U mx err
---|---|---|---|---|---|---
0.572063 | 3.409225e-05 | 1.023427e+00 | 1.023432e+00 | 2.404311e-03 | 1.034107e+00 | 1.034111e+00
0.286032 | 1.037969e-05 | 1.023429e+00 | 1.023430e+00 | 1.349427e-03 | 1.034108e+00 | 1.034110e+00

Order = 9

h | P comp | P err | Pmx err | U comp | U err | U mx err
---|---|---|---|---|---|---
0.572063 | 2.328387e-05 | 1.023428e+00 | 1.023431e+00 | 2.033288e-03 | 1.034107e+00 | 1.034110e+00
0.286032 | 7.124397e-06 | 1.023429e+00 | 1.023430e+00 | 1.141177e-03 | 1.034109e+00 | 1.034110e+00

Order = 10

h | P comp | P err | Pmx err | U comp | U err | U mx err
---|---|---|---|---|---|---
0.572063 | 1.664200e-05 | 1.023429e+00 | 1.023431e+00 | 1.746755e-03 | 1.034108e+00 | 1.034110e+00
0.286032 | 5.085321e-06 | 1.023429e+00 | 1.023430e+00 | 9.803705e-04 | 1.034109e+00 | 1.034109e+00